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ABSTRACT 


The problem of computing the temperature and heat flux 
distribution within combustion chambers is introduced’ 
and briefly discussed, highlighting the Hottel-Cohen | 
Zone Method. aa eer of localised errors in using this 
method are indicated with reference to systems having 
little or no flow recirculation. A suggestion is ad- 
vanced for minimising these errors by incorporating in 
the system model an equilibrium gas composition pattern 
that is explicitly linked with the jet-mixing process 
and the temperature pattern. The steps PSRaeHGeES the 
solution of the problem for the resulting complex model 
are outlined. Pursuant to this three computer programs 
are developed and shown to execute satisfactorily. The 
first program, EQUICALX, computes the equilibrium com- 
positions of the combustion product gases of any fuel at 
any combination of temperature and pressure. The second, 
GRAYGAS, computes a mixed multiple-grayt+tl-clear gas 
idealisation of the radiative behaviour of the burned 
gases of any hydrocarbon fuel. The third, FIFIELDX, 
computes the temperature and heat flux distribution for 
an axi-symmetric system using the output of the previous 


two for two operating conditions with gaseous propane, 
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C3He, fuel. A uniform gas composition pattern is pro- 
visionally adopted to facilitate the development of this 


program. The interfacing of these programs to facilitate 


the solution of the refined, complex model introduced is 


recommended for future research. 
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property of the real gas (see Chapter IV). 
total moles of component k at the start 
OP ereaction'. 

coefficients of pd in the polynomial 


expressions for the a 1/K, and 1/K¢ 


Gril 
(see Chapter V and Appendix D). 
concentration; mass iiraction. 
hypothetical concentration, mass fraction 
(respectively) in a thorough pre-mixture of 
primary and secondary flows. 
(c-C')/(C,-C"). 

Ay A 

a variable (see Appendix A). 

a mean constant-pressure specific heat 
(mass basis) - see Chapter V. 

the Craya-Curtet number. 

eMse e Wis hr aloe 

diameter ratio. 

diameter. 

emissive power (=0T"). 


a vector (equation (46)). 
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F formula for chemical species. 

aye al emissive power functions 
(see Chapter V and Appendix D). 

G, ea Gibby se function -sithe Sor tiaterot at 
Gibb* s; function: dG/on,. 

Gy. G/RT, G;/RT. 

G5, gs direct-exchange areas (gas-gas and 
gas-surface). 

GG, GS total-exchange areas (gas-gas and 
gas-surface). 

eG, Gs. directed-flux areas. 


h convective and bulk enthalpy flow heat 


tEransier coelreicients. 


Laie ssciila tee radiative equivalents of h (see section 5.4). 

H enthalpy. : 

sO O ; 

Tek numbers of species and components, 
respectively. 

Ge SG absorption coefficients, mean absorption 


cosetracvents. 


Lp li length, beam (path) length; L/54, 
L's L dimensionless quantity defined following 
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eciiatwOnaGc jee 

m total number of gray gases in an idealised 
multiple-gray-gastl-clear-gas mixture. 
(Chapter IV). 
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molecular weight. 

matrices. 

numbers of moles. 

total number of zones; number of surface 
zones (Chapter V). 

total pressure; partial pressure. 

energy flux. 

energy flux density. 

Badial"caistancevor a point £rom burner axis. 
universal gas constant. 

a vector (equation (46)). 

Reynolds number. 

surface-surface direct-exchange, total-—- 
exchange and directed-flux areas. 

a characteristic reaction time 
(non-dimensional) for the fuel-air 
combinacton. "y (Chapter «lL L)r 

absolute temperature. 

VELOGLtiy.. 
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Uunierstep function. (Chapter LV). 
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volume; volume flow rate. 

see Appendix A. 

distance s(axi1al) sofia pointe iromerebunmier 
(jet) mouth, also a variable; x/5.5- 


see Appendix A. 
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symbolyfor a general zone (G,S and g,s). 


absorptivity. 

the 'original' and the modified Thring-Newby 
Peci eulation Criteria; eresbectively. 

the Kronecker delta. 

a neighbour selector operator —- distinguishes 
between the neighbouring and the non- 
neighbouring zones for each zone, i. 
emissivity; change in value of specified 
variables from iteration to iteration, 
as defined following equation (23). 
wavelength. 

empirical constants in equations. 


(1) and (2). 


a step-size factor in equation (24c). 
chemical potential (see equation (12)); 
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UW/ RS 

cocEfievents (Chapter LLL). 

air/fuel ratio; theoretical air expressed 
as the. simple ratio of actual to stoichio=— 
metric air/fuel ratios. Alternatively, 
the inverse of the equivalence ratio. 


an error parameter, defined in equation (28). 
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mobecularn (as or dtretusivity = section 2.15). 
maximum. 

products (as for he in equation (3)). 
Rosseland-mean; reactants (as for oa in 
equation) (3)) 

at radial distance r- inethe redial. direction. 
related to the Rosseland-mean (Chapter V 
and Appendix D). 

Surtace; Surrounding fluid; \surtace—-ambirenc. 
stoichiometric. 

turbulent. 

at axial Gistence x; inthe axial direction. 
spectral (monochromatic) value. 

at the zero value of the coordinate in 
question (e.g., x, ©); Specified values 

0! Py - equation (19)). 
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CHAPTER I 
INTRODUCTION ; 


The design of combustion equipment often requires that 
attention be paid in some detail to the various aspects 

of the eomenea oven re ces S: These include the requisite 
size of the combustion chamber, the aerodynamics of the 
interior of the chamber as well as the transport of 
thermal energy between the hot burned gases and the cooler 
walls of the chamber. Considerations of the size of the 
combustion chamber, the combustion efficiency and thermal 
efficiency point to the necessity for a study of the aero- 
dynamics of flow within the chamber. Certain flow condi- 
tions, (for example, swirl) amprove the,supply of fresh 
oxygen to support combustion. This helps to reduce the 
size of the region occupied by the flame thereby af- 


fording some saving in combustion chamber size. 


Most combustion equipment consist, in essence, of a 
coaxal circular jet system in which the fuel issues from 
the inner jet. The annular flow comprises the secondary 
air supply. For non-swirl flows the proportion of the 
chamber occupied by the flame depends upon the nature of 
eae aa surrounding the flame jet. This may or may not 


be recirculatory. For non-combustion (i.e., cold) systems 
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some Similarity criteria have been developed to enable the 
prediction, of the existence or otherwise of. recirculation 

eddies as well as their location. For combustion, systems 

this task, as: much more difficult.. The aise: of the flows 

in these systems has naturally taken the findings from 


Ehose, Of cold syvsrems as ithe, starting, point. 


Heat transfer between the gases and the chamber walls 
depends upon the composition of the gases and on the tem- 
peratures of the gases and the walls. This sometimes 
necessitates some inquiry into the computation of the 
chemical equilibria of the combustion processes under 
various conditions of temperature and pressure. There is 
general agreement on the predominance of thermal radiation 
as the mechanism of heat transfer in combustion chambers. 
This follows from the fairly high temperatures usually 
encountered in such chambers and a realisation of the 
radiative behaviour of the burned gases, chiefly carbon 
dioxide (CO,) and water vapour (HO). Convection plays 
but a minor role while conduction, through the gaseous 


medium, hardly deserves even a mention. 


The passage Of radiation, through a. volume of gas,is)es- 
sential ly an. exponentials function of local. temperature, 
Bressure (the Combinedspartial pressure Of the “Tradvatively 


active! gases) —. Via) thesabsorptiiony coefficient (K) - and 
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path (beam) length. This is stated concisely in the fami- 
liem beer's Taw. Consequently, the transport of radiation 
through the non-isothermal, non-homogeneous systems usually 
encountered in combustion practice is quite a complex 
Subject.) JAS Of necessity, certain simplifying (often over= 
simplifying) assumptions have had to be made. The more 
immediate need has always been to have some estimates of 
the overall gas and wall temperatures and wall heat flux 


to enable questions relating to overall thermal efficien- 


cies, economics and metallurgy to be dealt with. 


The earlier analyses, therefore, have had to rely on some 
rather gross idealisations of the combustion systems. They 
also had to take account of the relatively primitive faci- 
lities available for carrying out the tedious calculations 
involved. Additionally, the dearth of thermal “alone 
data had to be put up with. Two of the most well-known of 
these simplified combustion system models are the Long 


Pornacesandethe Well-stirred Furnace [19,23]*. 


The Long Furnace model was usually applied to cases in 
Waien sche lateral. dimension, cf, the, furnace could be .con— 
sidered small compared with the longitudinal dimension. 


In such cases the gas and surface temperatures and the 


* Square brackets delineate numbered references detailed 
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wall heat flux could be considered as varying uni- 
dimensionally - i.e., in the longitudinal direction only. 
The combustion is idealised as essentially instantaneous 
and adiabatic. Use is made of the emissivity charts 
available for carbon dioxide (CO.) and water vapour (H,0) 
as well as the empirical relationships existing between 
the emissivity and absorptivity of these gases. The com- 
putation and the matching of the temperature and heat flux 


wiievds' proceeds alliong a trial and error Toute. 


The Well-stirred Furnace model was designed to fit such 
systems as justified the assumption bE completely homo- 
geneous and isothermal gas volume and a single wall. (sink) 
temperature. The implication of this is that the turbulent 
mixing taking place within the chamber was very vigorous 

and isotropic. Energy balance and some physical realism 
were allowed for by introducing a more or less stepwise 

drop in gas temperature at the exit. The information avail- 
able from this model was in the form of gross gas and wall 


temperatures, heat transfer and thermal efficiency. 


With the advent of the fast digital computer came the ever- 
Widening Opportunity to obtain more detailed information, 
atthe design stage, about a combustion system (be ita 
steam Plant boiler, an 01) retinery furnace, an oven, a 


laboratory-scale combustion test furnace etc.) than could 
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be furnished by the afore-mentioned models. More sophisti- 
cated models could then be constructed to better represent 
Ricerca l- iretsystem than. was Lormerly possible. — One of 
the more successful developments in this Ponce tion has 


been the Hottel-Cohen Zone Method [19,29]. 


In this model the gas volume and the confining walls are 
represented as assemblages of simpler, smaller-sized units 
called zones. Each zone is assumed to possess more or less 
uniform attributes within itself in the form of temperature 
(and hence black emissive power) emissivity and absorp- 
EIVILyY.. FOr example, a circular cylindrical, chamber may 

be divided up by a series of coaxal (with the chamber 
itself) circular cylinders and another series of planes 
perpendicular’to the chamber axis. This produces gas zones 
that are either whole cylinders in themselves (i.e., those 
on the chamber axis) or annular cylinders (farther from 

the axis). The surface zones may be short cylinders (on 
the curved surface) or circular planes or yet annuli (on 


the end walls). 


These zones are regarded as effecting the overall thermal 
energy transfer of the original system through simultaneous 
and successive interchange among themselves. The radiative 
component of such interchange is theorised as taking place 


through ‘exchange areas‘ existing between all possible zone 
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pairs. These exchange areas involve, in the first in- 
stance, the relative geometrical locations of each zone-pair. 
There are, on the whole, three types of exchange areas used 
in the analysis of radiative transfer - direct-exchange, 


total-exchange and directed-flux areas. 


The direct-exchange areas (ss, gS, gg) are, in essence, a 
measure of the proportionality of the net direct radiative 
exchange between any zone-pair to the difference between 
their black emissive powers. “The zone-pair is taken in 
isolation from the rest of the system. The value of this 
exchange area depends on the nature of the intervening 
medium, as well as the zones themselves if either or both 
of them are gas zones, the shapes of the zones and their 
relative geometrical locations. Tabulated values of these 
areas exist for cylindrical systems [19]. These tables 
use aS a parameter the product of the mean absorption co- 
efficient of the interval between the zone centres (K) 


and a typical zone dimension (B). 


The total-exchange areas (SS, GS, GG) additionally take 
into account the contribution, to the net radiative trans- 
fer occurring between any given zone-pair, of multiple re- 
flections taking place throughout the chamber. They are 


calculated from the direct-exchange areas. 
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Eachwor éthe Gdirected-flux areas (ace cop Gc) measure the 
proportionality which the one-way radiative flux passing 
from one zone to another bears to the black emissive power 
of the emitting zone. When the (real) aoe the inter- 
zone space is idealised as a mixture of gray gases, the 
directed-flux area becomes merely a weighted sum of the 


total-exchange areas corresponding to the different gray 


gas components. 


All the exchange areas are then seen to depend on the com- 
position and temperature of the intervening gas medium. 
This dependence would greatly complicate the application 
of the Zone Method if some simplifying assumptions were 
not possible. In combustion problems it is customary to 
ascribe a uniform composition (usually in the .completely- 
burned state) to all the gas zones alike. This implies 
thesuriformity of tthe gray gas absorption coefficients 


TREGUGHOUGRENe System, 


Using the appropriate exchange areas, energy balance equa- 
tions are written for each zone. These equations are non- 
linear inasmuch as they involve convective heat transfer 
between surface zones and their neighbouring gas zones, 
bulk flow terms for enthalpy transport between contiguous 
gas zones besides the radiation terms. The first two are 


linear functions of zone-temperatures while the radiation 
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terms are dependent on the zone emissive powers (i.e., the 
fourth power of the zone temperatures) and have temperature- 


dependent coefficients. 


Linearisation is sometimes possible whereby one can make 
somewhat easier the task of solving the equations which 
solution must necessarily proceed along an iterative route. 
The results of such solutions are the zone black emissive 
powers (E) from which the zone temperatures are easily 
derived. The distribution of the (surface) zone heat 


fluxes is then only a step away. 


Hottel and Sarofim [19] offer a partial example of the ap- 
plication of the method to a furnace in the shape of a 
rectangular parallelopiped with square end walls. Siddall 
[29] presents a step-by-step résumé of the method for a 
general furnace setting out a possible iterative procedure 
for solving the nonlinear energy balance equations. He 
does not give a specific numerical example of its success- 


tile appLication. 


The shortcomings of the Zone Method are the assigning of 
the Same composition to all. the, gas. zones of -the system 
and the assumption, of complete combustion. iene thonmed, peat 
will be recalled, is based upon the expectation that the 


turbulent recirculatory mixing occurring in the chamber is 
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vigorous enough. This expectation is not always satis- 
fied, as in combustion chambers operating without swirl 
engmunderssuch low Conditions that flow recirculation 

is minimal. Under such flow conditions the role played 

by molecular and turbulent diffusion phenomena in effecting 
the mixing process. is not insignificant. The very nature 
of these phenomena connote the existence of concentration 


(partial pressure) gradients in the chamber. 


Hottel and Sarofim [19] comment on the applicability of the 
uniform concentration assumption to situations where the 
Variations of the local concentrations of the radiatively 
active gases and/or their relative proportions are negli- 
gible. For other cases cited above the need arises for 
corrections to be made in the direct-exchange areas com- 
puted on the basis of uniform concentrations. Reference 
[19] suggests some apparently empirical corrections in- 


volving local absorption coefficients. 


Pieri et al. [27] discuss in some detail some of the im- 
plications of the uniform-concentration assumption. These 
concern the zone temperature and heat flux distributions. 
The discussion centres around a furnace in the shape of a 
rectangular parallelopiped of square section. They illus- 
trate some of the improvements to be expected under certain 


flow and operating conditions by allowing for the existence 
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of partial pressure gradients. These improvements vary 
in significance depending on the combustion chamber size. 
The assumption of a uniform concentration leads to the 
over-—estimation of the heat fluxes of the wall zones in 
the middle of the’ chamber as well as the corner wall zone 
temperatures. It also under-estimates the temperatures 
of most of the wall zones, the gas zones along the burner 
axis, the heat fluxes on the burner-end wall zones and on 
the exit-region wall zones. The errors in the temperatures 
can be as high as 600 F. The heat fluxes apparently 
exhibit a lower sensitivity than the zone temperatures to 
errors in the radiative properties of the system [22,27]. 
The discrepancy between the values computed with uniform 
and non-uniform concentrations can sometimes be as high 
as 110% of those calculated with the latter assumption 


[27]. 


The above-mentioned errors are attributable to the com- 
promising effect introduced by the uniform-concentration 
assumption between the high-concentration (of the combus- 
tion products) and the lower-concentration zones. Examples 
of the former category of zones would include the exit- 
region zones, those in the burner region typifying the 
Jatterecategory-melnids COmDromisingseriectlaleovexi sts 


between the higher-temperature zones (e.g., those in the 
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middie Of che furnace) and the cooler ones (e.g., those 

in the burner region). Where the concentrations are over-— 
estimated the (gas) zone emissivities will be too high. 

The zone temperatures are therefore more likely to be under- 


estimated [27]. 


What may be considered an alternative and no less realis- 
tic assumption postulates the existence in each gas zone 

of chemical equilibrium appropriate to both the zone tem- 
DeracUreroaid tie mole fraction, Of eLucl (on, alternatively, 
an air/fuel ratio) that would occur as a result of the jet- 
mixing process. This automatically introduces a pattern 

of concentration gradients coupled to the temperature field 
and the jet-mixing process. It would also account, im- 
Dplicitly, for fany reducing effect that dissociation would 
have on the zone compositions, and hence on one gas zone 
radiative properties. For lower-temperature zones the 


radiative properties of the gas volumes would be well ac- 


counted ‘for. 


Consequently, it would seem ideal to integrate with the 
Zone Method a zone concentration pattern in accordance with 
this postulate. This would seem to connote the updating, 
during the iterative solution of the energy balance equa- 
tions, of the zone concentration pattern. A refined, com- 


plex model of this nature should provide a more realistic 
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picture of the radiative interaction among the zones of 
the system, thereby improving the accuracy of the final 
computed distributions of temperature and heat flux. 
Hottel and Sarofim [19] suggest that Fils Beers labour in- 
volved in such a refinement of the combustion system model 
May not always be worthwhile. However, the rather large 
discrepancies cited earlier as existing between the pre- 
dicted temperature and heat flux distributions for the 
uniform-concentration model and those for the more refined 


variety suggest that such extra labour may sometimes prove 


eit aa 


The study reported herein is a preliminary step in that 
direction. Three computer programs are developed. The 
first Of these, —BOUICALX, isecapable Of scomputing the 
equilibrium states of reacting mixtures of fuel and air 
at various temperatures and fuel/air ratios. The pos- 
Sibility cL varying the overall syStemepressure sis —aiso 
allowed for. The second program, GRAYGAS, can compute a 
mixed multiple-graytl-clear-gas approximation for the com- 
bustion product gases of any hydrocarbon fuel-air mixture. 
The, last program, FH TEIELDX, predicts .the temperature and 
heat flux distributions within an axi-symmetric test com- 
bustion chamber using gaseous fuel. Use is made of the 
Hottel-Cohen Zone Method. The preliminary nature of this 


study dictated the provisional adoption of the uniform- 
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concentration assumption with the expectation that analysis 
OLethe more, Sophisticated model will follow from a:mating 
and extension of the procedures embodied in the computer 
programs and presented herein. Two operating conditions 

of the combustion system are considered and the associated 
flow conditions are. characterised by the appropriate values 
Of the Craya-Curtet Number (Ct). One of these is de= 
liberately chosen to be identical with that for the 

(mildly recirculatory) flow surrounding a free-jet flame 
burning in a stagnant atmosphere [5]. The other is much 
smaller, representing a more strongly recirculatory flow 


Situation. 


Hhewpsesentation covers the varios) aspects oF they combus— 
tion phenomenon pertinent to the study. These are separated 
into a number of more or less self-contained chapters, with 
the Appendices serving aS a repository for the requisite 
back-up information. The subject of combustion aero- 
dynamics is broached in Chapter Il, focussing on where it 
pertains to flame length and the use of flow recirculation 
parameters. Chapter III follows with an inexhaustive 
treatment of the theory of ‘complex chemical, equilibria in 
the gas phase. Next to be discussed, in Chapter IV, is 
therorocedure for Obtaining an idealised mul tiple-gray-— 
plus-clear-gas representation of a real gas mass. Chap- 


ter V details the steps undertaken to set up and solve the 
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applicable system of simultaneous energy balance equa- 
tions. Chapter VI offers some appraisal of the study 
undertaken and the results obtained, as well as some 
recommendations considered pertinent to the quest for 
improved accuracy in mathematical modelling of combus- 
tion systems. The Appendices, besides providing back- 
ground and eee aie information, also contain the 
listings of the three computer programs - EQUICALX, 


GRAYGAS and FTFIELDX. 
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CHAPTER® iL 


COMBUSTION AERODYNAMICS ' 


2.1 Flame Length 


The physical scat of a flame, whether pre-mixed or 
diffusion, is an important item in the design of combustion 
chambers. This is due, at least in part, to the un- 
desirability of the soot deposition that would inevitably. 
result from the impingement of a flame on the relatively 
cool walls of the enclosure - the wall quenching effect. 
Such contact between the flame and the chamber walls would, 
as can be imagined, have a deleterious effect on combustion 
erricwency.. As noted in Chapter £2, certain ~low condi— 
tions, swirl, for example, help reduce the size (length 

and lateral spread) of the flame, thereby permitting some 
reduction in the overall size of the combustion chamber. 
The optimum size chosen for the chamber must necessarily, 
therefore, take into consideration the expected maximum 


flame length. 


An early theoretical study of enclosed laminar diffusion 
flamesewas cneported 1nel928eby sBurke and «Schumann [7)- 
Numerous other theoretical and experimental studies have 
since been reported on free and enclosed laminar [3,4, 
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Naturally many different viewpoints have been represented, 


some more successful than others. 


For laminar flames the outcome of these investigations is 


that the length of the flame is a function of the volume 


aga 


content of that fluid. The latter variable can be repre- 


flow rate, of the burner fluid and the primary air 


sented by the non-dimensionalised air/fuel ratio, os (theo- 


retical air; ,,/ where Be is the stoichiometric 


ie 


alr tuel@ratio) Seeihesivame Pengti is, tosall “intents “and 


Sst. 


purposes, independent of burner diamter. Of the numerous 
expressions put forward relating the laminar flame length 
to the flow conditions and fuel properties, two stand out 


as the most widely successful. 


One of these is a semi-empirical correlation due to Hottel 


and co-workers [17,18] which is expressible in the form: 


— 2 * its 
ie = hs 10949 (Ve ath) + AG Ci 


in which ih L is the dimensionless length of the laminar 
f 
Uh uw 
ls A‘ and Ay represent- 
flamey Lb. Ge, (Le 1/94): the symbols , an 1 p 
ing empirical constants. The quantity tt Lsca ChalacteriSteac 
dimensionless reaction time for the fuel. It essentially 


measures the time taken, since issuance from the burner 


mouth, by the fluid on the burner axis to attain the 
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stoichiometric proportions deemed to exist everywhere on 
chew Elame envelope, ince luding the* flame tip | 7747 ..el. 
This characteristic time measure is — found. & -be 

{4 ln [(1+907 0) / (1408) 1}74 which shows its he dandedcd on 


o,- hts to ‘be noted that tt is a minimum when o5=0. 


This implies that the longest laminar flames are the 'pure' 


adifttusion fLlames. 


The other important expression for the laminar flame length 
fou cacOrecrcal one proposed Dy Wonlvev aie Pat) ee This 


Cay be pub an the forn: 
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ris 2 Oe rests ) (2) 


where Ay and ry! are constants whose values depend on the 


fuel and the primary air content of the burner fluid supply. 


The quantity Li L represents the dimensionless group 
J 


(Ve pte/ (a, Du ps Equation (1) is said to work best for 


flames characterised by high values of Ve B and whose 
’ 


lengths fall into the medium-large category. Equation (2) 


is supposed to be quite accurate in predicting the value 


Of he r over most of the range of variation of the latter. 
Ld 


Fauations (1) and (2))indicate the dependence of the laminar 


ie on V and oe ASoy. increases so 


Fou! £,B £,B 


does Reps the Reynolds Number of the (cold) burner-tube 


flame length, 
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flow. At a stage this enters the turbulent regime. The 
flame remains laminar, though, on account of the vastly 
InCEcased Gas Viscosity to be found inethe vicinity of the 
(hot) flame envelope. Beér and Chigier [6] report that 
increases of the order of 15 to 20 times the cold value 


are common. 


Naturally such an increase in gas viscosity can only be 
expected to have a considerable ‘laminarising' influence 
on the turbulent fluid flow issuing from the burner 
Orifice. At some critical value of Re, (turbulent), how- 
ever, the flame itself breaks down into a turbulent, 
brushlike affair - the turbulent flame. This breakdown 

is initiated at the flame tip and progressively approaches 
the burner mouth as Re, isincreased @1t se unable co 
reach the burner mouth, though, on account of the existence 
of a stable laminar stub at the base of the turbulent 
flame. This stub is referred to as the breakpoint length. 
Rabie, Lantvom reference 16) Plwsts tthe cri tical: Re, £ornea 
number of different fuels. 
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RAELE LT. ‘Critical Reynolds numbers, 


Reow(CoOld)., for wvartois. nels 
EES | 
Gaseous fuel Rex10 > 
Ho 2 
* = 
H, 525-6.5 
City. gas 3-4 
Gieyegaas Des Oe Oo 
(exe) 5 
Propane, C3H, S20 
Acetylene, CoH, , 
Methane, CH, 3 


* with primary air : 


In the turbulent-flame regime the relationship between the 
eS mer, we mc and the flow variables as well as the 

, 
fuels ipropéerties is) very di fierent (from tthat ain} theilaminar 


regime. y inetact, I. is independent of V for a given 


F,T £,B" 
burner diameter it is essentially constant. This constant 
value may be lower or higher than the maximum laminar flame 
length, depending on the) burner diameter. This is evi- 
denced in the experimental results reported in reference 


IiBa)\2e For smaller d., it is lower, in which case the flame 


passes through a maximum as Re, varies from the laminar to 
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ties cUrbulent tegimes. For Larger d,, the flame length 
increases monotonically even as the critical Re, is 
eraversed. All this is depicted in Pigures#la andvlb. 


As is the case with laminar flames, turbulent flames have 
Called rorth quite a few different viewpoints to their 
analyses. These approaches have led to a number of 
mathematical expressions relating the flame length to the 
various fuel properties and flow variables. Wohl et al., 
[34] extended their analysis of the laminar flame to the 
turbulent. They emphasised the diffusion phenomenon in- 
volved, taking into account its necessarily turbulent 
nature, and obtained a somewhat reliable formula. The 
usefulness of this formula is hampered by the need for 
arbitrary correction factors to be used in order to ac- 
count for the under-estimation of the effect of the eddy 


dat fusiv aby 


Yagi and Saji [35] based their analysis upon an even more 
detailed knowledge of the theory of turbulence. Using a 
guasi-cylindrical flame model and enunciating a ‘unit 
flame' concept, they obtained an expression for Ler 


which, though more complex than that of Wohl et al., is 


Of more Limited applicability. 
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Figure l b. 


Free-jet Flame 
Lengths vs. Gas 
Velocity at Burner 
Mouth. From 
eference [34]. 
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Hoectel@et ale 4,171) using the most general aoproach, and 
arguing from the point of view of the dynamic similarity 
of non-isothermal coaxial jets, vis A vis the relationship 
between the momentum change and buoyancy effects, have 
deduced, for situations characterised by low relative 
importance of buoyancy, a very general formula for Le op 
This formula has been shown to give satisfactory accuracy, 


averaging about 90% of the measured values. It can be put 


into the following form: 


T Ae Gus 
= PAS 70 BOSS 7/2 
L - x, = 10.6 Fol gay aced ' (1 tam! (3) 
ate PE 3 (Np Mp) st,cc Me ar 
in which Ch = 1/C,;," , being evaluated at essentially 


ScOlchniIOMmecr ic. conditions. ©The ratio np/np is determined 


for complete combustion. The quantity Xbp stands: ifor the 


dimensionless breakpoint length. This quantity was shown 


‘ 


to; bear A linmageneral, valimaximumeratio @fcabouc@0 419 tomthe 
Length tom) the saturbulentsportion ok) tie teuxbulent yil amey 


isa? (ep - Xp) ° This is true for most fuels, the ex- 
ception being CO for which the ratio has a maximum in the 
neuchhbourhoodsot.0, 84 el li ple lig this era tio egenemadsl yide— 


creases with increasing Re,- For most fuels the average 


Value tsi ‘about }0. 05. 
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A comparison of the formulas proposed by Hottel et al., 
Yagi and Saji, and Wohl et al., is presented graphically 
in Figure 2 taken from reference [17]. The universality 
of equation (3) is evident, M./M, being a parameter as 
would be expected on the grounds that gas diffusivities 


are strongly influenced by their molecular weights. 


Beér and Chigier [6] quote a formula for Le T that is at- 
7 


tributed to Guenther; 
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Leg = ie) (ome + 1) (De p/Pp) (4) 


the value of Pr being more or less the same for most gaseous 
fuels. This formula is alleged to be quite accurate (to 


Within 10% of thermeasured values). 


Enclosed flames in general behave differently from free 
flames. This is due to the effect of the confining walls 
on the flame aerodynamics. The confined flame is more 
prone than the free flame to instability resulting from 
recirculation. The confining walls set an upper limit to 
the amount of oxygen available to the combustion process. 
When the secondary air flow is too low, 1.¢6., when the 
WPeea ground stream’ [5] falls Short of the entrainment 
Capacity’ of the flame jet, flow instability is the likely 


result. Such instability usually manifests itself in the 
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Figure 2. A Comparison of Various Expressions 
for Turbulent tame Lengtns. 


From Rererence | 71. 
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form of flow recirculation which can therefore be thought 
of as an attempt by the flame to entrain the necessary 


fresh oxygen to support its combustion processes [5,33]. 


The confining walls also effectively prescribe a physical 
limit to the, lateral expansion of the flame. Confinement 
also increases the flame length besides tending to reduce 
the relative importance of buoyancy effects. Beér and 
Chigier [6] illustrate the increase in flame length due to 
confinement. The reduction of buoyancy effects is more 


Significant for horizontal than for vertical flames. 


Duwi li miow Recirculation 


Burke and Schumann [7] applied the term 'over-ventilated' 
to flames whose secondary air supply is more than is neces- 
sary to ensure complete combustion. Such a flame was 

shown to possess a closed envelope whose shape can be a 
prohibitively complex mathematical expression even for the 
very Simple cases treated in the quoted work. For a tur-— 
bulent flame this envelope exhibits a random oscillation 
about a mean position, being penetrated intermittently and 


randomly (in space and time) by turbulent eddies. 
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Wnen the secondary air flow 1s too low the Elame 2s said 

to be 'under-ventilated' and its envelope was shown by 
Burke and Schumann to fan outwards radially until it merged 
with the confining walls. This would indicate the com- 
plete entrainment, by the flame jet, of the secondary flow. 
This situation usually precedes the establishment, within 
therenclosure, of recirculatory flow. Figure 3 illustrates 


Ene» poine. 


Recirculatory flow is a complex mathematical problem even 
for the simplest isothermal, iso-density, free-jet systems. 
It has been possible to establish some applicable simi- 
larry tChicenia.!(|5.,0,32,53)..) [iis 1s. bDasedeupon theres— 
tablished fact of the reasonable constancy of the free- 


jet-spread angles. ; 


The treatment of the recirculation phenomenon for enclosed 
and free non-isothermal, non-isodensity, (combustion) 
systems has usually taken the form of an extrapolation of 
the results for the simpler isothermal, isodensity, free- 
Jatevarlety. | Sunavala etyal.§(32)_ point outs theyclose 
similarity between the behaviour of the early stages of 

the flame jet with that of a free jet.| (The conlinement or 
a flame jet necessitates the inclusion into the appropriate 
Similarity Criterion OL the. ratio Om thestypical lateral 


dimensions of the combustion chamber and burner, for 
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example, ‘the diaticter ratio, D, for a cylindrical system 


[33:5 


Arguing from a physical point of view, Thring and Newby 
[33) have shown recirculation to depend on the relative 
"strengths' of the primary (jet) and secondary (ground- 
stream) flows. This relativity is expressed in the ratio 


Xen! * 4p which exists between the axial distance, from the 
plane of the burner mouth, of the plane at which complete 
entrainment of the secondary air flow would occur (x) 

and that of the plane at which the steadily expanding jet 
would hit the walls of the enclosing chamber. This ratio 
can be shown to be proportional to the quantity CD, the 
first member of which stands for the (hypothetical) mass 
PLaction oL ther cued anya thorough pre-mixture of the 
primary and secondary flows. This product is essential for 
the definitions of the similarity parameters known as the 
Thring-Newby recirculation criterion (Bee and the Modi- 
fied Thring-Newby criterion (Bann) [6,32] (see Appendix A). 
Either parameter can be used to predict the occurrence or 
otherwise and the location of recirculation eddies. The 
'raw material’ forvall thas includes simply the mass flow 
rates of the primary and secondary flows and the dimensions 
Gertie wet.) (burner) and yconfining chamber {oj or com— 


bustion systems Bony has been found to be quite a reliable 


Simulacicy parameter. | This is depicted dn Pagures 4. inis 
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FIGURE 4. Use of Ban and Bony as 
Similarity Parameters 
for enclosed Jet Systems. 


From Reference [33]. 
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figure also illustrates that for isothermal systems Ben 


is the more appropriate flow parameter. 


An analogous parameter to 8, is the Craya-Curtet Number, 
Ct, which is derived from the fluid-dynamic equations of 
motion. It has been shown to be more or less identical 
with Bin within ees Of values from zero to unity 
[6]. Being derived from the fundamental equations of 
motion, it is the more rigorous and accurate indicator of 
flow similarity [5,6]. Becker et al. [5] have indicated 
the nse of ¢Ct. in eens isothermal, isodensity jet 
systems and have described it as a unique criterion, for 
a given diameter ratio, D, of dynamic similarity. It 
usually varies in value, for jet-type flows, from 


fp 


EOe gs tally DOR aL Ven Ds EOL recirculatory 
flows they have shown that the Ct can be used, contrary 
to the suggestion of Beér and Chigier [6], to statisti- 
cally correlate the different attributes of the recircu- 


lation eddy - its eye, upstream edge, downstream edge, and 


the aesociatea volumetric flux cf ftluid. 


Becker et al. also discussed the Significance of the Ct 
value of 0.75. It represents the borderline between flows 
that are recirculatory (with Ct smaller than 0.75) and 
those that are not (with Ct larger than 0.75). Flows whose 


Ct fall between the values of 0.5 and 0.75 possess very 
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jicvtele recizculeation. Smaller values of Ct represent flows 
having progressively stronger recirculation. For such 
flows the assumption of a uniform concentration is tenable. 
On the Other Nand, for flows that have lvttle- or no re— 
circulation (2.e., with Ct greater than about) 0.5). the 
validity is increasingly questionable. Becker et al. also 
demonstrated quite conclusively the possibility of re- 
producing in an enclosure flow conditions that are, to 

all intents and purposes, dynamically similar to those 
associated with a free jet issuing into a stagnant free 


stream. Such flow conditions are characterised by a Ct 


O£f4'0. 6:73. 


For any given system (with a given D and fluid supply con- 
ditions) there exists a 1:1 correspondence between Bon 

(ox Bam? and 6' (the theoretical air) as well as between 
Clhuandpos, (see Appendix, A)Layror the systemostudied,, 0 

head aevalue. of 72 which Gave values OL (Cl Of 0207/3 and 
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CHAPTER LiL 


CALCULATION OF EQUILIBRIUM COMPOSITIONS 


Gigli Vbevetalevekbrenmabe)s| 


The radiatively active gases of principal interest in com- 
bustion studies are carbon dioxide (CO,) and water vapour 
(H,0) - The emissivity as well as the absorptivity of any 
given volume of combustion product gases depends upon the 
partial pressures of these two gases, among other things 
[19]. The uniform-concentration assumption, throughout the 
combustion chamber, discussed earlier in connection with the 
Zone Method (Chapter I) has been portrayed as leading to 
localised errors in the calculation of temperature and heat 
flux distributions. These errors are more serious in the 


former tenan the Latter [22,27] « 


Gas zones can attain temperatures that are high enough for 
the effects of dissociation to become signiticant.. “inethat 
case the zone concentrations should be computed to take 
account of the phenomenon. Thus the radiative interactions 
of the zones of the system will be more truly assessed. An 
analogous situation could arise during, the siterative ssolu— 


tion of the energy balance equations mentioned earlier. 
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The analogy would become clear upon considering the itera- 


tive solution to be something of a non-steady-state operation 


of the system, the iterations representing some hypothetical 
EimessCaleo Thus the validity Of the tinal esults) of. the 


computation would be enhanced. 


An additional motivation for the investigation of the com- 
putation of the equilibrium compositions of reacting systems 
was the need to have on hand a documented facility of that 
nature. Its inclusion herein also lends completeness to 


Civs study. 


Combustion and chemical literature abound in studies of the 
Scuttibrla OL reacting systems (12,137 260, 50) , 2Zelezniiand 
Gordon [36] provide an extensive list of references. Their 
paper deals in some detail with the chemical, thermodynamic 
and mathematical aspects of the subject, and it is their 


methodology that was adopted and will be briefly outlined. 


In general, there are two schools of thought on the problem 
of “equilibrium computations: the’ S“equilabrium constant tore 
mulators' and the 'free energy minimisers" [36]. By es- 
tablishing the relationship between equilibrium constants 
and free energy, Zeleznik and Gordon have shown the two. 


parties to have much in common. Differences between them 
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seem to lie chiefly in the detailed manipulation of the 
equations describing the reacting systems. 

In essence, the computation of chemical equilibrium can be 
visualized as the solution of a set of non-linear equations 


of the form; 
Q'(x) = 0 (5) 


where x is vector representing the state variables of the 
reacting system. In the equilibrium constant approach these 
equations are essentially equilibrium constant expressions 
for elementary reactions involving the species present in 
the system. In the free energy minimisation approach the 
determination of the extremum of the free energy function, 


Q(x) say, becomes the SOLluUCiLOn Of the equation; 
(d/dx) 2 (x) = 0 (6) 


A comparison of equations (5) and (6) indicates two pos- 


sible relationships between 2'(x) and Q(x): 


~ 
ie 
iH} 


(d/ax) 2 (x) (7) 


Doe) = (Lye) ee Mae (8) 


Mr being any PD matrix, the unit matrix, fOr ianstance, and 


qr ft) the transpose of 2' [36]. 
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Naturally the solution of the equilibrium problem would be 

greatlysracilitated if the equations (5) are kept as simple 
as possible#,; This is achieved through the use of the iMor— 

mation reactions of the chemical species present in the 


System. 


3.2 An Overview of the Theory 


For any system operating at a constant temperature and pres- 
sure, the equilibrium state is marked by the fact that the 
Gibb's free energy of the system is at a minimum. This con- 
Sideration would apply to a large number of combustion 
systems. This free energy function is a function of the 
system temperature, pressure and composition. In other 
words, for a chemical system containing ing species the 


Gibb's function can be stated as 


Gi=G (T,P,n, Mor ++ Tio) (9) 


It is an extensive property and as such must satisfy the 


equation 
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The equilibrium condition can therefore be stated as, 


- 
Dome ail =) 0 ‘ial 


For gases, to which the theory presented here is inten-— 
tionally limited, the chemical potentials (us) are defined 


as 


= * 
wy us (T) tORT in (n;/n) P (72) 
where the us (T) are the zero-pressure chemical potentials 


at the temperature T and n is related to the n; by 


a Sen (13) 


In writing the formation reaction equations mentioned sais 
the reacting system is considered to be no more than an 
assembly of sa species of which ee essentially constitute 

the 'building blocks' for the remainder. These k° species 


are referred to as the ‘components’ of the system. 


In a 'C-H-O-N' (Carbon-Hydrogen-Oxygen-Nitrogen) combustion 
system the appropriate species could include C,H,,0,,N5,CO, 


CO. CNH HCO,HO.,H.0,O0,0H,;N;NO (1° = 15) “cf which the first 


Le 2k 


four constitute the components (k° 


I 


4). These are the 
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chosen sets of species and components for the combustion 
system considered in this study. 
The formation reaction equations for the remaining foe 


species are then of the form; 


= 0 Kee ie (14) 


where the Vane are stoichiometric coefficients and the = are 
the chemical formulae of the species. These formulae can be 


mathematically symbolized in terms of the components as, 


ko 
reas AT eG) i er tee <P 3 (15) 


iuewhach v.. us the number of molecules of the ae component 


ki 
contained in one molecule of the alae species. The coeffi- 
cients Vio and Vig must then satisfy the relation, 
{oO 
2; Viiv a = 0 (16) 
ul J 
Gumi: 
alc kO 7; k <jS1 
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The chemical law of conservation of the atoms (molecules) 


of the components has to be satisfied at all times during 


the chemical reaction. 


A conservation equation can therefore be written for the 


molecules of each component giving, 


2 (een P= b=" U7) 


where by is the number of molecules of component k at the 
start of the reaction. The use of equations (16) and (17) 
Ebanstorms therequilibrium) condition (equation) (11)) into 


the form; . 


Valk as ae 


the dimensionless form being preferred. 


To complete the thermodynamic description of the system it 
is requisite to specify any two state variables chosen from 


among T,P,S,H and V. In combustion problems the choice 
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might be, 


and Jae ea i (19) 


where Po and T> are constants, 


Bauartonse( lS) gm (id), ek ls) and, (19) constitute a system.of 
ee equations in as many variables (n being one of them). 
In theory a solution is possible. The system is non-linear 
on account of equations (18) so the solution must eroeees 


iteratively. This is considered in the next section. 


36050 Method, of }Computation 


A suitable procedure is the so-called Newton-Raphson itera- 
tion [36]. In this method some initial estimates of the 
variables x are made, and at each subsequent Lteration. £ 


the current values of x are usea to recalculate the variables 


for the next -- (I+l)th -- iteration in a functional manner, 
Pew; 
ie a 20 
Sse ED as ony, J (x,) 2 (x7) ( ) 
= a0" (21) 
or JT(x,)-AXp 14 Q" (xy) . 


Wite rer: is the Jacobian matrix and AXs44=Xy417%1 
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This Kindson functional ateration can: in theory. be carried 
COULSON ines system of, equations, (13) ,- (lv)e (18) anawti9)) 


thereby progressively approximating the equilibrium com- 


position of the reacting system. However, one important ob- 


servation needs to be made. In theory, as the equations 
stand, there is nothing to prevent the mole numbers be- 
coming negative at See stage in the iterative process. But 
physical realism and the very nature of equation (12) (which 
should be always satisfied by the variables in equation 
(18)) demand that the mole numbers be always seaieives An 
effective artifice Chae eines this prospectively awk- 
ward situation is to use the logarithms of these quantities 


as the variables rather than the quantities themselves. 
Using this approach and performing one or two convenient 


variable transformations (see Appendix B), one would obtain 


the following set of equations: 


Bee a) tb Alnene ct vps G, (22) 
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ee 
2 (E.t.) + (=e) Aln n = ei + c (23) 
k'=1 
Leak cake | 
40 
t — 
where ie = 5 (Vy iene ne) 
i=l 
Ups = Aln Nee + ete Ain 7 
.O 
a 
= ] 
ek e Vea 
i=1 
,0 
ask =e ee 
iSplis k oe Kae: 
i=1 : 
© 
a 
and En ee yaa 8; n; 


In an expanded form these equations look like this: 


1 i ' { 
Pl), aie a © 440 MS, Uy E byes l= 
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oY 2S 2 Us 2H eae 4 
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EXO, §0, +--+ 0,0 £&,0 u,,0 Ey 40 + SURSHIEN 
aT bo eu E,9 (-e_) Aln n —e +G 
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At each iteration (I) this matrix equation is easily solved. 
The mole numbers and the total number of moles are updated 


uSing the U; and the approximate relation: 


ore 


ce 
fj 
3 
»* 
R 


“X1) /X, (24a) 


from which Xray * Xz + x,Ain xy (24b) 


The risk of over-correction, which would be run should equa- 


tion, (245) be used as 1t stands, is considerably reduced by 
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using a modified version of the same equation; 


X74] 2 X, + ) a x Alin Xr (24c) 
in which A''' is an empirical step-size factor. Reference 
[36] .recommends eee ae be determined such that, 

JAln n] < 2 (25) 
Aln (n, /n) < 2; for (n,/n) = ey © yee) 
Aln (n; /n) eee bce (n,/n); fox (n,/n) < ‘iy = (27a) 


Equation (27a) may be written more informatively in the form; 


8 


Bm =in (17a); for (n/n) <10— (27b) 


Aln (n,/n) < dr (On 


Evidently the value of A‘''' varies from iteration to itera- 
tion. Its use, subject to the conditions listed in equa- 
Grons (25) ,. (260) and (2/a) or (275), would, limit the growth 
of all currently significant species (with n,/n 2 me =) to 

a maximum of aye (equation) (26) in. ine CuULTenELysdnscign iii 
cant species (with n,/n < way Ss are also prevented from sud- 


denly rising to mole Fractions in excess of Ta (equation 


(2a) OL Z 1) jie 
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Tie i eeratiVer process 15 Carried on until acceptable Limits 
Of error “are attained cfiihe? approach toi those slamits: vs 


monitored by an error parameter (wW) given by, 


12 Ae 
li Z : 
Wee Lee) 4 ef pee y/? (28) 
Gee shegh 
; a8 eee 
where e— =u 2 Vi... It can be seen from equation 
Hy 1 k,i'k 


(28) that » is essentially something of a 'root mean square 
error', following each iteration, in the computed species 
mole numbers, total number of moles and the chemical poten- 
tials. When this quantity is small enough, depending on the 
desired accuracy, the problem of computing the equilibrium 
composition of the reaction products will have been expended 


Woh 


3.4 Computed Results 


Appendix B contains a listing of EQUICALX, a computer program 
written in FORTRAN IV to carry out the computation procedure 
described above. Sample results of that program are re- 
produced in Table II for heptane (CH, ¢) atee /00GR tor ichres 
values of @'. Similar results for the same fuel are taken 


from reference [8] and included in that table for comparison. 
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TABLE If. Mole fractions computed by EQUICALX for the 
products of the combustion of yeptane 
(CH, ¢) iNvaLe ac 27007R- (1 500ek) stor 
Clase), Sore Ormand 200 

compared with 

Similar results given in reference [8]. 


Computed by EQUICALX from Reference [8] 


0.00001 
0.00085 
OS S516 
0.00002 
0.12408 
050 


0.00000 


0.14181 
0.00000 
0.00002 
0.0 


0.00008 


0.00000 
OeLOd 
0.76006 
0.00000 
0.06436 
0.0 

0.00000 
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0.00000 
0. 07/352 
0.00000 
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0.00089 


03 24559 


0.00000 
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The effectiveness and accuracy of the above procedure are 
now evident. Discrepancies between the concentrations of 
Ny in corresponding columns of that table almost exactly 


accounted for by the difference between pure and atmos- 


pheric No. 


TABLESE UL S eMokeatractiions kof iGO ,<CGO. rand so 0, present sin 


2 2 : 
the products of the combustion of C3H, with 
air in the temperature range 1000 - 4500 R 


KO OD 
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TABLE IV. “Mole fractions moieco ; CO. and H,O present in 


the products of the combustion o£ C3H, with 


air in the temperature range from 1000 R to 
4500 R (d¢"' = 3.0) 


HO 0.0545] 0.0545; 0.0545] 0.0545 |0.0544 /0.0539 0.0524 |0.0484 


It is possible to tabulate, using EQUICALX, the gas composi- 
tion following the combustion of any fuel at any temperature 
and theoretical air (¢'). This is done for gaseous propane 
(C,H,) in the temperature range 1000 R to 4500 R and at two 

Values=om 0). (1.0 and 3.0), as cdeplcredgin Tables ttl and 

3 and H50 are listed 


Since the latter two are of the most interest in the study 


TV Only the mole. tractions Om CO; CO 


reported. CO is included merely for comparison with the 


other two gases. 


The integration of the equilibrium computation scheme into 
the calculation of the distributions of temperature and 
Walleheatetilux is a faimly strargnttorward matter, as wild 


be discussed in Chapters V and VI. 
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CHAPTER: EV. 


THE RADIATIVE PROPERTIES ‘OF COMBUSTION PRODUCT GASES 


451i Introduction 


Thermal radiation is a form of energy transport whose emis- 
sion or absorption by matter is a consequence of the pos- 
session or results in the increased possession, by matter, 
of temperature. It collectively refers to the infra-red, 
the visible and the ultra-violet forms of radiation. 
Thermal radiation represents, therefore, a small portion 

of the much wider family of electromagnetic wave phenomena 
whose members range in wave length, A, from the very-short- 
Wavelenguhay says. Chrougn the X-rays, sthermals radiation, co 
the long-wavelength Hertzian (radio and electric) waves. 
These wave phenomena have at least one thing in common: 
they all travel at the same characteristic speed, the speed 
of light. Consequently their frequencies vary in a reverse 


manner to their wavelengths as the spectrum is traversed. 


Electromagnetic waves result from or give rise to energy- 
level. transitions. withinethe atomic/moleculazn structure of 
matter. The atom is theorised as consisting of a central 
core, the nucleus, surrounded by a number of electrons all 


statistically distributed in discrete energy levels, the 
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number of electrons per nucleus depending on the species 

of matter (atomic number) and its state of ionisation. 

The nucleus itself is composed of even smaller units than 
itself, including the so-called eon akira rode and 
neutrons) and other more esoteric particles. These nucleons 
are also arranged in nuclear energy levels. Atoms bond to- 
gether to form molecules which may be symmetric, e.g., 

2" Ho, SO, anadsCo. 


These molecules are perpetually in vibratory and, in fluids, 


H510, and No, Or FaSymmMetric, e.g., CO 


rotatory as well as translatory motion. The vibratory and 
rotatory motion of the molecules take place 


only in a number of discrete associated energy levels. 


It is the transitions that occur among the energy levels 
enumerated above that are responsible for, or ,that do 
result from, the interplay between radiation and matter. 
These transitions occur only with the emission or absorp- 
tion of radiation possessing energy exactly equal in magni- 
tude to the energy differentials among the energy levels 


in question. 


The Quantum Theory postulates electromagnetic wave pheno- 
mena to be lumps or particles (quanta) of energy. Each 
quantum is associated with a particular wave (frequency) 
and with a particular energy-level transition within 


matter. The y-rays are associated with nuclear energy-level 
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transitions, X-rays with those of the lower electronic © 
energy levels, whereas the ultra-violet and the visible 
parts of the electromagnetic spectrum are involved in the 
transitions occurring among the higher electronic energy 
levels. The vibrational and rotational ener crat eval 
transitions occur sequel to the emission or absorption 

of infra-red Bee ee The longer-wavelength members of 
the spectrum, e.g., the Hertzian waves, have to do with 
the excitation of (i.e., energy level transitions among) 
the valence electrons of matter. Such excitation is the 
result of the application to matter of an electric poten- 


tial. 


Thermal radiation is emitted and absorbed by all bodies 

at all: temperatures: above! absolutewzero. ~The amount emitted 
and absorbed by matter increases with the absolute tempera- 
ture of the source (emitter). Generally, in most computa- 
tions of engineering heat transfer, thermal radiation be- 
comes important only when the absolute temperature of the 
emitter is thousandsvo£ degrees. Thermal radiation), as 
hinted to above, is dispersed through the spectrum of wave- 
lengthsieethisyis duehtorthetstatnsticaltenaturesor all 
atomic movements. The general level of such movements de- 
termines the general "level' of wavelengths of the signi- 
ficant emitted radiation, as depicted in the spectral- 


distribution (spectro-radiometric) curves sketched in 
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Pegure 2 fOr an ideal (black body) Tadtator. ) The ‘black-— 
body radiator emits (and absorbs) the maximum possible 
radiation at any given temperature. The shapes of these 
curves are given by the so-called Planck's Law [19,24]; 


CoAT 


B, = (C)/A7)/(e ay) (29a) 


r 


where Cy and C5 are the Planck's constants (first and 


second, respectively) and E, the monochromatic emissive 


r 


power of the emitter. 


ies is evident from Figure 5 that, in general, the hotter 

a body, the shorter the wavelength range for the major 

part of the emitted radiation and the more peaky the spect- 
roradiometric curve. The converse is true for cooler 
bodies. For example, the sun, at a temperature of about 

10 500 R emits over 90% of its thermal radiation between 
0.1 and 3u whereas a body at about 2 500 R emits between 
Peand 201 (24.7 The totalythermal radiation emittedeby a 
body is the area under its spectroradiometric curve. This 


has been found to be related to its temperature as follows; 


ee teil! | (29b) 


which is the Ste £an-Boltzmann Law, o being the Stefan- 


Boltzmann constant. 
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For solids, the vast majority of which are opaque to thermal 
radiation, this phenomenon is essentially an affair of the 
solid surface. Most solid surfaces do not emit thermal 
PadlaviOn hn CONnLOLrmLiey With tie Reh ee ae curve 
of the ideal radiator. The monochromatic emissive power 
of conductors at a given temperature decreases with in- 
creasing wavelengths whereas it increases irregularly for 
non-conductors. The monochromatic emissivities of the 
(real) surfaces vary in a similar fashion with wavelength. 
It is usually possible, however, to adequately approximate 
a real surface with a 'gray' surface - one with an emis- 
SiVicy peat is Constant througnouc thesspectrum 8197.24), 
The spectroradiometric curves of gray surfaces then dupli- 


cate that of the black-body surface on reduced emissive- 


power scales. ‘ 


All the statements made above as well as those to follow 
Willestali hold at fom ememissivity’ One read) “absorpci— 
vity', bearing in mind that the temperature of interest 
is, and will be, that of the emitter. Use of the emitter 
temperature ensures the satisfaction of Kirchhoff's Law 
which requires that for a body in thermal equilibrium 

, = €)- For a gray body this requirement becomes a = &€ 


r 
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Unlike what happens in solids, thermal radiation in fluids 
is a volume phenomenon. In combustion studies, gases are 
of much greater interest than liquids, from the point of 
vView,of. radiation heat..transfer..«The emission vef+ thermal 
radiation in gases (as well as its absorption) is a 
‘function’ of the nature and the composition of the gases 
as well as the temperatures of the emitter and the absor- 
ber. The emission from the gas proper is non-luminous, 
but the presence of solid particles like soot adds a lumi- 
nous component to the radiated energy besides introducing 


the phenomenon of radiation scatter. 


The Bouguer-Lambert Law states that the attenuation of a 
collimated beam of radiation passing through a volume of 
gas is an exponential decay function. This leads to the 
following expression for the monochromatic emissivity of 


the gas 


(30a) 


where the Ky is the monochromatic extinction coefficient 
for absorption. ».,.Theassunption, es, tacut hy madesbhat, the 
gas is homogeneous, in which case K, is constant along the 
beam length... This beam-length, L,, is usually defined) as 


the radius of a hemisphere containing the gas whose radia- 


tion, Ls, werceivedual the «centre of pthejilaiwihace oF In 
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general, K) varies throughout the spectrum. An ideali- 

sation is usually conceived in the form of the 'gray' gas - 
one whose absorption coefficient, by analogy with the emis-— 
Sivity of a gray surface, is constant throughout the spec- 


trum. It is a very useful concept as will be seen shortly. 


A real gas has neither a spectrum-wise constant absorp- 
tion Coesficient,, nor One that is a/continuous: function, of 
wavelength. In combustion the term ‘real gas' usually 
denotes either co. or H,0 or some mixture of the two. 
These gases have radiative properties that are complex 
functions of the partial pressures of the gases, the total 
pressure and temperature. They are subject to several 
"broadening' phenomena including pressure-, collision-, 
pelr— andi linée[proadening. | They also voverlap in several 
wavebands, resulting in their radiative properties not 
being strictly additive. In addition they selectively 
emit and absorb in some wavebands whilst being transparent 
to the rest of the spectrum. This selectiveness has 
origins in the resonance of the molecular rotations and 
vibrations that have been referred to earlier. A conse- 
quence of this behaviour is that a plot of K, versus 1/2 
usually is a series of isolated roughly triangular shapes 
dispersed through the spectrum, ac sketched in Figure 6. 
The values of Ky are more or less symmetrical about the 


mean values of 1/A (corresponding to the apices of the 
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triangles) The gaps between successive 'triangles' are 

referred to as the 'windows' and they collectively repre- 
Sent, tne gclear Foas (radiatively transparent) component 
of the real gas. Charts have been compiled and published 


for the emissivities of C05, H,O and other gases. Hottel 


2 
and Sarofim discuss, these at some length in their book. 
Some complex correction factors and schemes have been de- 


veloped to reduce the errors introduced by the various 


phenomena cited herein. 


It is evident then that some simplifications are mandatory 
if the emissivity charts of the gases CO. and HO are to 

be usable for the purpose of solving the energy balance 
equations cited in Chapter I. Fortunately this is possible 
in the form of an idealised multiple gray gas ‘and clear gas 


mixture that approximates the radiative behaviour of the 


real gas. 


4.2 Gray Gas Approximations 


The basis for this artifice can be established quite 
simply by taking the definition of the total hemispherical 


emissivity of any volume of gas. This is, 


= E, da 
eg F Jean Eyed he nN 
O O 


7 ~ a ; ; er) 
‘ ; ae a _e x a . : 


2 


feston hia een iii secinvnales > re 

dase ausiy ck sispnet smote iets eauciekh atrial Ga >: 
-sh nsed ved domod>a bus 226563 nokissz10> xelgmoo emom 
avoiisy edd yd bsoySotsni exerts. aft sovbex of beqoley : 


. 
yaotebasm ous anoliasoitiiqmia smoe Jars neds jnebive ei aI n 
od 815 Og Ens 09 esesp sit 30 atzedo ytiviestne odd Bt 
sonsisd Yproms sit palyfoe to Ssaeqzng sft 202 eidseu od 

eldteacg et aint ylotenvsx10i -.1 reagend nt hetio enoiseups 

esp t69f5 boo esp yore Sigqisinn bseiissbi as to msot sdt ak | 
eid 20 qWoiveded svidsibs: sii astemixorgge tedd exutxim 
«esp {sez 


znoitemixorgygA end-ystd Seb 


etiup borleifdstes od neo ecitisys aint yo2 siesd edt 
Seotzeiiqeimet (stot edd te noitintieb ons pnides yd baci 
vai ekdT asp 30 sawkov yas 26 yJiviseime 


* *) VE Pa “er 
oO” Oo 


a 


. awe 


which can easily be transformed into 


1 
ua 2 eed UE, 
eae = (1 e lads (30b) 


r ; 
where f = | [eny Pec. This quantityeic tne ew raccion of 
the total emissive power occurring within the wavelength 
interval from zero to A. It becomes unity, naturally, when 


che Uppers Limi trOoL integration, 1s) taken ‘to Infinity. 


By invoking the congrvence between an integral and the summa- 
tion of discrete strips of area under a curve, - 
it is possible to re-express the integral in equa- 


tion.(320b)ias thessummation: 


) (30c) 


in which the K. are constant, the Af. being the fraction 
of the total emissive power falling within the range(s) 
of X} to which the mean of K; + AK/2 would apply (see 
Figure 6). The number (m) of the discrete items under 
the summation cee can conceivably be very large. It is 
evident from Figure 6 that the ranges of X to which 

K, + NK72. apply are not neccessarily contiguous. However, 


this is really immaterial as long as the proper walues of 
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Af. are used. Using the definition of a gray gas, equation 
(30c) can be looked upon as the weighted sum of the emis-— 
Sivities of a number (m) of gray gases. The weighting 


FaccLors, Aft; Can be alternatively symbolised asa 


Chine 
that equation (30c) becomes, 
a3 
= Syeaan 
ca cee (l-e ) aes 


It has been found in practice that m rarely exceeds 3 [19, 
20,22,27]. For this study m = 2 was found to yield satis- 
factory accuracy. The temperature variation of aan demands 
that the fan and the K be functions of temperature also. 
The potential difficulty resulting from this is obviated 

by keeping the K. constant at some mean values while the 
aed are forced to account for the temperature variation 


Of € The procedure for computing the a and the K, 


g,L’ G,i 
have been outlined by Hottel and Sarofim [19] and will be 


briefly reviewed. 


4.3 Method of Computation 


Mesemissivity Of any irregular gas volume, as) percerved at 
a point within the volume, is a function of beam length, 
Lew this varies With cirection. Lt is customary, thererore, 


to define a mean value, La This is recommended to be 
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Using this mean value, a table is prepared of EG 1 Versus 
: , 


L while the latter varies from the minimum value to be en- 
countered in the system to something more than twice the 
mean beam length [19,29]. This is done at the temperature 
in the middle of the expected temperature range. Equation 
(30d) is then recast in the more general LOrm: 

i~l m 
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where U' (i-1) is a unit-step function having a value of 


@er@ (Ores! and unity Lon all 4-1 in ovner words, =1t 


is zero during the computation of ac 4 and Kye 
a 


(the gas emissivity at infinite beam length) is 


The quan- 
EL y te 
ag 


foe) 
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given by, 


m 
€ Si -Ade eo (33a) 


where aa o 1s) the wei giernoetactom 1Or, Lhe rc learsgas Con 
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ponent. 


Tf the RHS of Gaquation (32)) is assumed tojbe arranged such 


that 
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: ie 
Ken =e Kerly forvawk bx (33b) 


Clem one fings that at “large *values Ot & tne RHS Of that 


t 


equation will be dominated by its first member, i.e., by 


=—Kei 3 : ie ; : 
ag ; & a". Consequently if a judicious choice is made of 
4 


EG | then a Semi-log plot, using common logarithms, of the 
? ; 


LHS of equation (32) versus L would, at large enough values 
Ob i, be a straight line of slope “kK, with an intercept on 


tie Vertical axis of logsa.”.. 
G,1 


A repeated application of this procedure will yield all 


the pairs of a and K, deemed necessary to satisfactorily 


Gri 


approximate the emissivity of the real gas at the given 
temperature. As indicated above the number of gray gases 
required is generally less than 3. The weighting factor, 


aa or for the 'clear gas‘ component is calculated residually 
’ 


from equation (33a). At other temperatures in the expected 


range the a are re-calculated (keeping the K; unchanged) 


Cy 2k 


by matching the computed with the actual values of emis- 
sivity at conveniently chosen values of L. For m=2 the 


suggested values of L are Lo and 2Le 


Having calculated the a at the selected temperatures 


‘Gypak 


covering the expected range, the next step is to obtain 
the best polynomial (ordinary) an 2 tO8iitl he temperature 


VaAclobdon OL Gach. “a A good way of doing this 2s first 
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CO TObcain fhe best Chebyshev polynomial (with the minimax 
DroOvercLYy fof such polynomials) tae fit each eye over the 
temperature range. These Chebyshev polynomials are then 
easily converted into their equivalent ordinary polynomials 
thereby retaining the accuracy of the former. Standard 


subroutines for these polynomial operations are available 


in the IBM SSP (Scientific Subroutine Package). 


The ordinary polynomials representing the ag 4 can be 
LA 


generally written thus; 


Jax 


2 enh ae Her 34a 
eae Pa (34a) 


where the ae are the Lemperature coetticuents, Denne being 
the degree of the polynomial. Siddall [29] suggests that 

in general a quadratic eae! would be adequate. However, 
for the cases studied (propane gas-air combustion at 

6' = 1.0 and 3.0) it was found necessary to use polynomials 
with 5) Snag to obtain an acceptable accuracy. This was par- 
ticularly true of the very last gray gas component to be 
computed which would suggest significant round-off errors 
accriling, £rom, tne Computation. Ol. the Ssarnlier eordy, Gas 
Components. The dotted portions of the, curves, for sek in 


Figure 7b and a in, Bogure So Llustrare This lasts poutine 


G,2 
Equation (34a) can now be used to express the emissivity 
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(or absorptivity) of the real gas mass, at any beam length 
(L) and temperature (T) as, 
5 ' 
m max KL 


€ (or a yaa slat ee Toe ¢ 
g g eae ee ) (l-e ) (34b) 


A computer program named GRAYGAS was written to carry out 
the above procedure and 1s listed in Appendix C. Sit in= 


corporates calls to the SSP subroutines mentioned above. 


4.4 Results 


Figures 7a (curve 1) and 7b present the output of this 
program for a test case. They show a mixed 3-graytl-clear 
idealisation of an equi-molal mixture of CO. and HO at 
2500 R. The total pressure of the gas mixture is one at- 
mosphere and the mean beam length of the enclosure is 

3.4 feet. The range of beam lengths used was as in Figure 
6-26 06 trererence [19], the tull (4367, the non—-—dotted) 
curve Of which as reproduced ian: curve ss On Figure wa. 

vs. L as computed by GRAYGAS inter- 


g,L 
polating within the HO and CO. emissivity tables. The 


Curve 2 shows the ¢€ 


latter were constructed using the emissivity charts in the 
Same reference. The discrepancies observed among these 


curves Can be atlLtrloutea to, a number of factors. These 
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include the extrapolation of the (semi-log) charts to higher 
temperatures; the interpolation within the small-scale 
Charts, especially at the very long and the very short beam 
lengths; the uncertainties inherent in une ceonrect tons efor 
the pressure-broadening and spectral overlap phenomena as 


well as the round-off errors inevitable in single-precision 


erithnmetic. 


It may be noted at this point that the procedure as set out 
by Hottel and Sarofim [19] requires m to be specified at 
the start of the computation, to be changed if the accuracy 
obtained is not acceptable. This need not be the case if 
it is borne in mind the assumption of equation (33b) in- 
dicates at least an order of magnitude difference between 
successive gray gas absorption coefficients. This connotes 
an order of magnitude difference between the weighting fac- 
tors. Furthermore, the maximum contribution of a gray gas 
to the emissivity of the mixture is its corresponding 
Wwelgntingatactor. Consequently,, the, neca to, calculate the 
ith (say) gray gas would not exist if its maximum contribu- 


ti0nLef.0 67 estimated sat. about is negligible. 


i 
Casal 10 Poea—14 
This ensures that no more gray gases are calculated than 
are necessary to approximate the real gas system to an 


acceptable degree of accuracy (See Appendix C). 
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The program, 


GRAYGAS, was used to approximate the products 


OL propane-gas-—air combustion at ®' = 1.0 and 3.0. ‘The 


results are shown in Figure 8 and Table V. 
lustrates the variation with temperature of the a 


the latter lists the hs 


TABLE V. 


and the Kx. 
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The former il- 


G,i while 


Gray gas extinction coefficients (K,) and the 


temperature coefficients for polynomial (in T) 


approximations to the gray gas emissivity 


weighting factors within the temperature range 
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CHAPTER V 
TEMPERATURE AND HEAT FLUX DISTRIBUTIONS - THE ZONE METHOD 
5.1 The Combustion System Model 


The system under study is now oes aine described. A 
horizontal cylindrical test combustion chamber 6 feet 
long and 1.5 feet in diameter encloses an axially posi- 
tioned turbulent flame. The flame is maintained by a 
constant supply of gaseous propane (CHe) through a 1/4- 
inch diameter fuel jet surrounded by a uniform annular 
secondary air flow. The mass flow rate of the fuel was 
fixed at 9.75 lb/hr. and two different air supply rates 
were specified so that $" = 1.0 and 3.0. The ambient and 


fluid supply temperatures were taken to be 537°R. 


This system is modelled as a closed cylindrical chamber 

with a plug flow directed towards the ‘exit' end and ex- 
hausting through an exit hole 6 inches in diameter located 
Oiietnesaxas. [ft 1S reduced to a collection of Zones, das 
and surface, in accordance with the Hottel-Cohen Zone 

Method. The gas volume is divided into 4 contiguous cylin- 
drical zones coaxal with the chamber. The curved surface 

of the chamber 1S Similarly divided into 4 shorter cylinders. 
The inlet end is taken as the fifth surface zone while the 


exit end, excluding the exit hole, makes up the sixth. 
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This Simple coarse-—zone configuration, sketched in Figure 9, 
was adopted in order to facilitate the development of the 
computer program required for the exercise discussed in 
this chapter. Naturally ,\a Liner zoning eaves be used and 
would be required for more realistic and accurate results. 
Further comment is reserved ror later this section. For 
each zone an energy balance equation is written. This 


equation can, assuming steady-state conditions, be summarily 


stated in the general form; 
rate of energy gain - rate of energy loss = 0 Gicy, 


The search for the solution of the simultaneous equations 
of this kind constitutes the mainstay of the Zone Method. 
For a typical gas zone, the first member of equation (35) 
would involve the rate of energy generation within the 
zone. Thus the combustion (heat generation) pattern needs 
be known if the problem is to be solved. For the axi- 
symmetric system under study, the chamber-fuel-jet diameter 
TAtIOfUS 72. Le seems Weasonable, when, io presume ‘the 
combustion process to occur along the chamber axis. In 
other words, if the system possessed more than one gas 
zone in the radial direction, then the energy generation 
term in equation (35) would not exist for all gas zones 


except those situated on the axis of the chamber and 
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within the length of the same occupied by the flame. 


The length of the flame is thus also required. It may be 
computed using one of the two relationships Leese eka in 
Chapter II. For this study equation (3) was used on ac- 
count of its established universality and reliability 
a EY als sEecieees a flame length of 4.7 feet. Equa- 
tion (4) was employed to check that computed flame length 


with which it was found to agree within 5%. 


The computation of the radiative exchanges among the zones 
requires information on the gas zone compositions. In 
Chapter I it was stated that in the application of the Zone 
Method a uniform, constant gas composition was usually 
assumeditec: bei G@ifect P19 727429) . “*Thirs usually leads 

to localised errors in the computed distributions of tem- 
perature and heat flux and these errors are sometimes 
sizeable [27]. Some improvements would naturally result 

if allowance is made for the ee sen pressure gradients 


occurring in the chamber [27]. 


It was suggested earlier that a zone-temperature and 
mixing-process-dependent equilibrium composition may exist 
in each of the gas zones. This supposition would strengthen 
somewhat the coupling among the gas zone radiative proper- 


ties, the zone temperatures and the mixing process. During 
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an iterative procedure the zone temperatures would vary 
from iteration to, iteration. The assumption of a zone- 
temperature-and-mixing-process—dependent equilibrium com- 
position would, therefore, add both to the ‘rigour of, and 
the computational effort required by, the combustion 


system model. 


It is often rewarding to imbue a theoretical model such 

as this with all available conceptual and mathematical 
PAgOUr eet the daight of .thisi, Hottel and Sarofim [19] 
have suggested the initiation of a theoretical modelling 
exercise with the most rigorous approach as a possible 
means of abstracting the value of the model, presumably 
through the subsequent, progressive and selective shedding 
of various aspects of that model. One cannot cavil at 
this proposition. However, it sometimes happens that other 
considerations make the opposite approach more attractive 
and feasible. Such is the case with the work reported 


herein. 


The setting up and successful Operation of the subordi— 
nate, constituent aspects of the complex theoretical model 
was regarded as a necessary initial step towards their even- 
tual merger into the composite system. Consequently the 
uniform, constant composition assumption was provisionally 


adopted with the expectation that the outcome of this work 
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will be of value in the eventual inclusion, into the Zone 
Method of analysis, of variable, non-uniform gas composi- 


ions. 


The program GRAYGAS (Chapter IV) is used to obtain an 
idealised multiple-gray+l-clear-gas mixture for the burned 
gases. This aoe also fits a polynomial in T to each 
Sains within the temperature range 1'000 R - 4 500 R as 
outlined. The output of this program serves as data for 


the solution of the energy balance equations to be dis- 


cussed later in this Chapter. 


5.2 Exchange Areas 


The real significance of the gas zone compositions lies in 
their effect on the radiative properties of the gas zones 
and, through the latter, on the exchange areas construed to 
exist between all the zone-pairs of the system [27]. In 
the non-isothermal systems usually encountered in combus- 
tion, the directed-flux areas are of predominant interest 
(see Chapter I). They are calculable from the total- 


exchange areas as follows, 


(36) 
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where Z represents the symbol «for a general zone, gas -(G) 
optsurface a(S) ; i, being unity for all gas-surface, 
surface-gas and gas-gas interchange but zero (to include 
the"clearegas tcontrabution)tior surface-surface inter- 


Change. “The welohting ‘factors ta are computed as in 


Gy. 
Chapter IV. The numerals i and j refer to the zones. The 
total-exchange areas (ZZ) are computed from the direct- 
exchange areas (zz) [19]. The latter have been computed 


for a number of combustion system configurations for which 


formulas and charts have been developed. 


Fore’eylindrical’chambers fHottel SandeSarofime [E9}lfiurnish 

a set of tables attributed to H. Erkku. These tables are 
directly applicable to chambers whose lengths and diameters 
aremsuchtthatia, commons typicaiezeneudimensiong_Basays ris 
easily stated. The tables then provide the inter-zone 
direct-exchange areas for various values of KB, where K is 
a gray gas absorption coefficient applicable to the system. 
There are as many direct-exchange areas, for any zone-pair, 
as there are gray gases in the medium between the zone 
centres. The tables are applicable to systems with up to 
12 ‘axial@and 5S radialigas zones. ’Theysthereforé,enable 


the directed-flux areas to be computed. 


One drawback on the use of these tables is that (cylindri- 


cal) combustion chambers are not always of dimensions that 
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are integral multiples of a common dimension, B. Further- 
more, Systems are often encountered with a larger number 
of axial gas zones than are covered by the tables. This 
usually results from attempts to keep KB within the range 
Of thet tables, (OP to de25)%0 Hottelt and Sarofim [19] tsuggest 
some approximate, but not so simple-to-use, formulas for 
such systems. Another apparent defect in the tables is 
the need to apply to the tabulated exchange areas correc- 
tion factors whose values are often unknown or difficult 
to determine. Additionally, Hottel and Sarofim advise the 
avoidance of KB values in the range 0.4 to 3.0 which prac- 
tically includes the range of KB used in the tables. For 
KE in excess of 3.0 (i1.6., for “optically thick*® gases) 
they recommend the diffusion approximation for radiative 
transporc 11,7 10,19,21,729). The values of KB encountexed 
im his Werk ranged from 1. Sco 359. ca (Secerlaples Vi =ror 


the values of Ke au ne cLoererent oj, 


This concept of diffusion is attributed to S. Rosseland. 

Lt considers the transport of radiation through an opti-— 
cally thick medium as being akin to the diffusion process 
Or to, thermal conduction. The ditftusing on conducted 
material is the thermal radiation quantum, the photon. 
Using this concept, formulae are obtained for the directed- 
flux areas. lt 1S ustial to assume that Che zones are. only 


influenced by their neighbours except in surface-surface 
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exchange in which the clear gas contribution is undimin- 
ished. 

Wie use Of the datiusion approximation iseitacilitated by 
thes det initarentort al certain’ absorption coerricient,, the 


ROS Setandemean@absorprtion=coetiicient, KL, which is given 


R?’ 
by, 
tyre, ee 
(1/K») = | KeOhey! SP hae (37a) 
fe) 
which can be put in the form 
uy 
(1/K}) = (1/K,) aie (37b) 
fo) 


where K, is the monochromatic absorption coefficient at 


wavelength A, £* being defined by 


r 
£* = | 2,728), ar Carsy 


the quantity £* is easily Shown to be relaved (oO 1 detined 
in Chapter IV (see Appendix D). With the tacit assumption 
of a step-wise relationship between K and A [19] and hence 


between K and £*, equation (37b). is then transformed into 
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HT 
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on 
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(1/K}) T Sac, i)/K, | (37a) 


which shows that Kn can be calculated from the Ac and 
? 
K. and as such is as much a function of temperature as 


they aren 


For two adjacent zones i and j whose interface area is 


Pau the directed-flux area is 


Tei Wee ey. (KB) (38) 
Legs ps ge eget F 


where T, is the temperature of the emitting zone and B is 
the zone centre-centre distance. This formula is directly 
applicable to homogeneous gas-gas exchange. For gas- 
surface, surface-gas and dissimilar gas-gas situations 
(2) ae gives erroneous results because of the 

ex clusion, fron its derivation, of the phenomenon variously 
known as ‘energy jump' [10], ‘temperature jump" [1,19] and 
‘radiation slip’ [21]. This phenomenon manifests itself 
in a temperature discontinuity at the interface. The origi- 
nal concept of the diffusion approximation, which works so 
well in the interior of a homogeneous gas mass, is unable 


to cetect, or account for Chis discontinuity. 
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Deissler [10] discusses this phenomenon in some detail. 

He cites some earlier works that use a 'first order energy 
jump' giving insufficiently accurate estimates of the tem- 
perature jump. He also proposes formulae for his ‘second 
order energy jump‘. which yield what appear to be a remark- 
able improvement in the accuracy of the predicted tempera- 
ture discontinuity in the region around the interface, a 
solid wall, for example. For an emitting-absorbing gas 


flowing uniformly through a tube of diameter, d he puts 


17 


forward the following formula for the directed-flux area 


(per unit interface area) between the gas and the wall: 


2 
-g Ee] Kee 
= -1 - r,W —1 We eS ' L ak “ 1 27 R S 
g § Ly ions Lee ad 
(39) 


where ay ana E. are the emissive powers at the tube centre 
and at the wall surface, respectively, Ke is related to 


the Kh and is defined as 


Motel iene me 
(1/Ko) = J /85) (98,798) da (37e) 


which leads to 
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da 


- at 


G, i) /Ke (S74) 


| H 


m 
(1/Ke) = p> (ag 3+ 
1=1 


The development of this relationship is outlined in Ap- 


pendix D. 


Peis nNoOcable that Kp and Ke are computed at the tempera- 
ture of the emitter to account for the dependence of es 

on the emitter bemperature. This requirement also ensures 
that Kivenore ve law is satisfied should the emitting and 
the absorbing zones both attain the same temperature. 
Equation (39) can be applied in the Zone Method if qd, is 
replaced with Ar, the radial width of the gas zone adjacent 


to the wall, Be and Es becoming the gas and wall zone black 


emissive powers. 


The radiative exchange between a flowing gas and the end 
Wablston a-cylindrical@chamber, taking account of the 
energy jump, has received no treatment in the available 
literature. As a first approximation, however, the formula 
developed by Deissler [10] for a stationary gas mass be- 


tween parallel walls may be cautiously adopted, viz.; 


= wi eI 
(GS sa) t= (qH22) = ZF Z(Kh Sax) + (2 - 5) (40) 
Ss 
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WheEreRAX ~1Sathe*axial-width' of =the gas~zone adjacent to 


the end = zone; Ey and Es being the gas and end-wall zone 


emisSive powers. 


For optically thick systems, the st paeeee trace inter- 
change may be approximated by that which takes place 
through the 'windows' of the radiation spectrum, i.e., 

the tcleanqas’ contribution. From equation (26) this would 


be 


(41) 


5.3 The Energy Balance Equations ‘ 


For a surface zone, the first composite term of equation 

(35) is made up of 

(1) the direct one-way flux from all zones, including 
self-irradiation and 

(2)4the internal convectivesheats: transtexeixom theenéigh— 
bouring gas zone, whereas the second term includes 

(1) the radiation emitted by the surface, of the general 
form ecA, and 

(2) convective thermal losses from the exterior of the 


system to the surroundings, in the general form 
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hA(T;-T)), which contains the applicable heat transfer 
coefficient (h), zone surface area (A) and temperature 
(T,) and the ambient temperature (T,)- 

If the zones of the system are so numbered that the first 

N‘' zones are surface zones, the rest being gas zones, 

then the energy balance equation for a surface zone may 


be symbolized in the form 


N' N N 
—SE ——>= 
y Sesie. + 2 G.S.E. + = Cie ane Ty Dan) 
oe _— I _— | 
acti 2) I psy Die eals ep eet a) aa 


(42) 


aac Done Tectn\tietoben te 


Lonel i i< N’, where N' and N are the number of surface 


zones and the total number of zones, respectively; See 


is a neighbour-selector operator, being unity for gas zones 


adjacent to the surface zone but zero for all others, he 
and ne are the gas-side and ambient convective heat trans- 
fer coefficients, Ais the interfacial area between zones 


Peace |p. and As the area of zone i. 


For a gas zone the rate of energy gain, neglecting conduc- 


f10on through the gas, will include 
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(i) 


(2) 


(3) 


The 


(1) 


(2) 


(a) 


The 


direct one-way flux from all zones, self-absorption 
inclusive, 

energy generation rate on account of combustion pro- 
cesses taking place within the zone, and 

the bulk enthalpy inflow from the neighbouring gas 


zones. 


rate of energy loss will be made up of 

all radiation emitted by the gas zone, of the general 
form 4KV, where V is the volume of the zone, 

convective heat loss to the neighbouring surface zones, 
and 


bulk enthalpy outflow to the neighbouring gas zones. 


plug flow assumption stated earlier has the effect of 


forcing each gas zone to receive bulk inflow only from the 


gas 


zone immediately upstream of itself. Similar con- 


aTderatlons would hold for the bulk oneriow. 


Li the gas zone numbering is carried in such a manner that 


the zone numbers increase downstream and outwards radially, 


then the energy balance equation for a typical gas zone 


takes the mathematical form: 
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for Nit < i < N, where m* is the mass flow rate per unit 
area, C. the appropriate mean specific heat at constant 
pressure for the temperature range encountered and oe 
the rate of combustion-generated energy release within the 
gas volume. For this study ont was zero for all zones 
except those residing on the chamber axis and containing 
at least some part of the flame. It was assumed to be 
directly proportional to the proportion of the flame con- 
tained in the particular zone. Equations (42) and (43) 


then constitute a set of N simultaneous equations in as 


many unknowns, vis., the zone emissive powers. 


5.4 The Solution of the Equations 


These equations are nonlinear, containing terms in E with 
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temperaturesaecpendent coefficients as well as terms in’ T. 
Their solution must therefore be obtained iteratively. 

Some simplification may be attempted to facilitate the 
Solution. “This may be in the form of a linearisation of 
the equations = by converting the terms in TI to their rough 
equivalents in terms of E [19]. When this is done, one 


obtains the following set of equations; 
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where h' Sp dadyd vd ene 4 A fee tag Aa E.-E. 
S,J s 5/ 7 g | J i)7 J i? 


rt = 0 (AI = = 
and aes m Gee Ae E,) all of them computed at the 
current temperatures and emissive powers. Bee is the inter- 


face area between neighbouring zones j and i, while A; is 


the area of the ith surface zone. 
The set of equations can be summarised in matrix form thus 
ME=R (46) 
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a vector comprising the Fas and R the vector of the RHS 


of the set of equations. 


Siddall [29] suggests a procedure for solving these equa- 

tions. The essence of that procedure is adopted, with 

slight modifications, and will be briefly outlined as 

follows: 

(1) fill in the vector R using the specified ambient tem- 
perature, convective heat transfer coefficient (ho), 


and combustion pattern, 
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(2) guess the distribution of temperature among the zones 
of unspecified temperature with which to 

(3) calculate the directed-flux areas (ZZ) in accordance 
With eequations (38), (39), (40)- and (41). Using these 
and the specified gas-side convective heat transfer 
COCLTICIenNL, 

(4) compute the elements of M and 

(a) 2SOlve fon E, the new temperature distribution. This 
is done either directly by Gaussian elimination or by 
the Gauss-Seidel procedure depending on the conditioning 
of M which in turn depends on the ordering of the system 
of equations as well as the current temperature distri- 
DuULIOn. & baving computed the vector E, then 

(6) test its members for agreement or otherwise with the 
previously computed values of the same. If there is 
considerable error, go back to step 3 above and repeat. 
If not, then one would have obtained the required tem- 
penbature Gistribution., For this, study an ,gcneement 
within 10° R was considered adequate. 

(7) Calculate the net heat fluxes through the surface 
Fones — the heat flux distribution. As va check on the 
calculations, 


(8) comonte’ the total rate of heat loss to the surround= 


ings as, 
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(FA) 


and the net heat generation rate (Q,) as 


= foe os Wal Bi 
Q,.= 0 mh Oe eit Tinlet? . (48) 


Where Qa is the overall heat release rate, An the cross- 
sectional area of the chamber, the subscripts to the tem- 
peratures being self-explanatory. The quantities Q and Q. 
should agree within acceptable limits. For this study an 
agreement of 80% or better was considered adequate for the 
coarse-zone configuration used on account of the approxima- 


tions pertaining to the composition, flow and combustion 


patterns and the energy balance equations themselves. 


A computer program, FTFIELDX, was written to apply the pro- 
cedure reviewed above. This program is listed in Appendix 
D. It uses, as part of its input, the output of GRAYGAS 
(see Chapter IV) to facilitate the execution of step 3 of 


the procedure. 


5.5 Computed Results 


Tie output of the program, for the assemblage of 4 gas 
and 6 surface zones, is listed in Table VI and is diagram- 


Watically displayed in Figure 10 for 0> =) 1.0 and' 2.0. 
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The curves in Figure 10 exhibit trends in general agreement 
with those reported from experiments [33] with, and theo- 
retical, predictions, [27]. fom, some .enclosed combustion 


systems. Figures lla, b and c illustrate these results. 


TABLE VI. The temperature and heat flux distributions, 
the Qa and the Qr, of the combustion system 
(oly =S1% Offand?s .0) 


end zones 
2 


all zones| 3 
4 
5 
6 
gas zones if 
8 
Ss) LB6o ok 1494.7 
yRG LOO9.L 1148.0 
Q 


2 (10 sBtuvHr.) 
Ry, 


. eytete: 80.8 
Oo (per cent.) 


Temp. (>) R) Heat Flux 
(Btu/hren ft<4n OR) 


2345.0 
3054.6 


bs gligl his ae 2048.3 
Loess oS SR aT 3340.7 
1494.9 5747.4 4107.9 


O31. 0 2364.2 26652 


1000. 1 
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It is clear that the iterative procedure set out above and 
employed in the computer program, FTFIELDX, yields results 
in qualitative agreement with the expected system be- 
haviour in practice. The maximum temperatures compare 
with the adiabatic flame temperatures. From a detailed 


quantitative viewpoint the temperatures and heat fluxes 


may not be in very good agreement with measured values. 


This is partly due to the assumption of a uniform, constant 


concentration and partly to the apparent contradiction be- 
tween the plug flow assumption and the pheonomenon of flow 
recirculation. The former source of error would be more 
significant for the low recirculation case (Ct = 0.673, 

@' = 3.0) while the latter would affect the higher recir- 


culation case (Gte= 0.230). 0i@= W.0) toa greater degree. 
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Figure 10. Temperature and Heat Flux Dietri bution 
Computed by FTFIELDX for the Confined 
Turbulent Flame Jet System Studied 
(O° Yer 0 eand 330) 
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Figure lla. 


Measured Axial Temperatures 
for Two Bnclosed Turbulent 
Jet Flames. Reference [33] 


Figure 1b. 


Computed Axial Gas 
Temperature Distribution 
for an Enclosed Turbulent 
Jet Flame (Uniform Con- 
centration). Reference 
[272 Rum +e 


Figure ic. 


Computed Wall Heat 
Flux Distrtbution 
for an Enclosed 
Turbulent Jet 
Flame (Uniform 
Concentration). 
Reference [27], 
Run 274 
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CHAPTER Vi 
CONCLUSIONS AND RECOMMENDATIONS | 


The coaxal jet system, which typifies most combustion 
systems, can be succinctly represented by similarity para- 
meters of the Were or the Craya-Curtet variety. 
These parameters encapsulate information about the aero- 
dynamics of the combustion system interior, vis Aa vis 

flow recirculation and jet mixing. The presence of strong 
recirculation within the combustion chamber (i.e., at low 
values of these parameters) would validate the assumption 
of a uniform composition throughout the chamber. For 
other systems without a strong recirculatory flow some 
allowance needs be made for the existence of concentration 
gradients in the interests of accurate predictions of the 
temperature and wall heat flux distributions within the 


chamber. 


A suggestion has been put forward of a means of making 
such an allowance, namely, the existence of equilibrium 
gas compositions, in the gas zones, that is linked to the 
zone temperatures Sal the jet-mixing process. Pursuant 
to this a computer program, EQUICALX, has been developed 
to carry out the computation of equilibrium compositions. 


This program has been shown to yield quite.accurate results. 
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The computation of the temperature and heat flux distribu- 
tions extant in the combustion chamber is greatly facili- 
tated if the 'real' gas contained therein is idealized as 
a mixture of a relatively small number of pee gases and 

a 'clear' component. To this end a computer program named 
GRAYGAS was developed and, as indicated by the discussion 
in Chapter IV, gives results that are satisfactorily 
accurate. This program provides a vital link between the 


emissivity charts and the computation of the temperature 


and the heat flux distributions. 


The execution of this latter step was carried out, for the 


system studied, with the aid of a third computer program, 


FTFIELDX, whose output was shown to agree qualitatively with 


published measured values. As was indicated, ,improvements 
in the accuracy of these results are to be expected if the 
theoretical model is further refined. This is due to the 

coupling existing among the temperature and heat flux dis- 
tributions, “the concentration gradients .and «the jeét-mixing 


process. 


Hottel and Sarofim [19] have suggested that when the varia- 
tions in the local concentrations of the radiatively active 
gases and/or their relative proportions are significant, 
the theoretical model should allow for the non-uniformity 


of the gas concentrations. The assumption of an 
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equilibrium gas composition pattern linked to the zone 
temperatures and mixing process is one way of doing this. 
The output of EQUICALX shown in Tables Pee een ma Ly. 
indicate that although the equilibrium SE fora 
given %', does not change very much with temperature as 
long as the latter is not too high (i.e., less than about 
SUUURRetOr O° = 120 and about 3500 (RY for 'v=) 320) ee tne 
relative proportions of CO, and H,0 change significantly 
with @'. The connection between this fact and the equili- 
brium gas composition postulate mentioned earlier should 


now be obvious. 


The procedures outlined in Chapters III, IV and V which 

are embodied in EQUICALX, GRAYGAS and FTFIELDX should prove 
useful in setting up and analysing the more sophisticated 
model that would account for this interconnection among 

the temperature distribution, the heat flux distribution, 
the gas composition pattern and the jet-mixing process. 

Wiis may take form of first computing the sVvalveror so. for 
each gas zone in the absence of combustion and as deter- 
mined by the jet-mixing process. The one-to-one corres- 
pondence existing between ®' and the equilibrium gas com- 
position of each zone at any given temperature should deter- 
mine the radiative behaviour of each gas zone. By using the 
procedures described in Chapters IV and V, the analysis of 


the model can be carried through to obtain the temperature 


and heat flux distributions. 
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APPENDIX A 
THEORETICAL AIR AND RECIRCULATION CRITERIA 
A.1 Relationship between 6' and Ct 


Becker et al. [5] give the following expression for Ct, 


a eS Bags aL al 2 pera 7 6 
Ce Uys IU. ees 12 oe De omueed (A.1) 
= # 1 ; 
where UL. hahaa rO yeep: ry + Loko and D is the diameter 
ratio dp/d » by MakingO@the Substitutions ¢ =) ib —— 
B 2 
Ce 
2 = : 
B= 1 °— 1/D? and U = U5 ,0/95,0 , one obtains 
3 ee ee A) pee ee) (A.2) 
o Dp? 


This directly relates Ct and U for any given D. Beér and 
Chigier [6] give an expression for the ratio of the secon- 
dary and primary mass flow rates that can be put in the 


forms: 
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For a hydrocarbon fuel CyHy» 


e e M Y 
M /M. = 4576 Me (X + 7) ot (ACA) 


Consequently, equations (A.2), (A.3) and (A.4) furnish a 
relationship between Ct and 6'. This is graphically de- 
picted in Figure A.1 for two temperature ratios 

eto 3/0 = 1.0 and 1.4) of the secondary air to the fuel 


jet, D serving as a parameter. 


From this chart Ct and 6° can be easily correlated. For 


example, if D=72 and Te = 1.0 (for the system’ 


10! <0 
studied) then a Ct lof 0.673 corresponds with an o* of 


3.0, while $' = 1.0 would correspond with Ct = 0.230. 


. ° t 
A.2 Relationship between $6 and Bae Bon 


Beér and Chigier give a definition of the Thring-Newby 


criterion that can be expressed as; 
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Ct 


1x1o7! 
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at Various Values ofr De 


ime modified Thring=Newhby criterion is defined [6] as; 


He 0. 
B aa caer eee J,O L/i2 
TN oan (— ) (A. 6) 
<7ies Ps,0 ; 
Hence 
- =,l 172 
Pen = Boon HMOs 7 - 1)] (Ado 
18 


Equations (A.3) to (A.7) inclusive embody the relationship 
between Bn or Bony and O°. Similar considerations apply to 
this relationship as were discussed in connection with the 


Craya-Curtet Number (section A.1). 
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APPENDIX B 
TRANSFORMING THE EQUILIBRIUM COMPUTATION EQUATIONS 


Be Llothes Equations 


Bquations (13), (17), (18) constitute a set of ae + 1 


equations in as many variables: 


¥ ~° - 
eee eV =O (18) 
a ee | aa ER 
OS eo: s ae 
= @, 
= 0 | ey 
5 =} 
x (vi sn.) + by fe) ( 
1=1 
Ie keer ae 
and 0 
i 
Te th. = O (13) 
; a 
1=1 


As commented on in Chapter III, these equations are non- 
linear and as such can only be solved iteratively. A 
recommended line of approach is the Newton-Raphson method 
[36]. Furthermore, the risk of obtaining negative mole 
numbers necessitates the adoption of some artifice to 
ensure the avoidance of that potential difficulty. This 


is by way of replacing the mole numbers with their 
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In applying the Newton-Raphson functional iteration (equa- 


EO (2)))) 


following; 
—- _-, 
Uy weit in ns 


from which follow 


ou, /2 (in n;) 


and 


du,/d (1n n) 


One then obtains the 


colt dae gt 


n, (du; /en,) 
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of equations: 


to the set of equations, use would be made of the 
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4° 
FeV Ale a bevallen = ¢! (B. 1b) 
eat ae. al Dyk 
ees aes ae : 
-O 
aL . 
n Alnn - Ratniees ns ae (By ic) 


where the e's are as defined in Chapter III following equa- 
tion (23). Upon inspection, equation (B.la) would seem to 


indicate the need to define a new variable; 
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BYEialkiicethat SUupStILUCLOn, ;equadtiGns 22 )esand (25) 60L 
Chapter III are obtained, the solution of which constitutes 
the object of the computer program, EQUICALX, listed in 
this Appendix. At each iteration the corrected mole numbers 
are computed using the Alin ns in the format of equation 
(240) e(Chapter Tit). These vain n; are in turn computed 


from the u; thus; 
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B.2 The Computer Program, >) EQUICALX 


A listing of this program follows. The program, written 
in FORTRAN IV, was run on an IBM System 360/67 machine at 


the University of Alberta Computing Services installation. 


(a) Tay (a) TA) ©) 


OAA OO) OA @ @e OGG @ @& oe @ OO @ 


ar OO CGY Ce @ Miia @) dal Gar (2 


EQUICALX 


COMPUTES EQUIL COMPOSITIONS IN COMBUSTIGN PROCESSES 
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METHOD OF COMPUTATION = NEWTON-RAPHS3IN FUNCTIONAL 


ITERATION AS DEVELCPED BY 
ZELEZNIK AND GSRDON 


DIMENSION STDMU(30+20)+CHEMU(S30) > 

« I MPC20) SSPECS( 30) sVNU( 3055) +QKZ05), 

; BNJ (3053) »sBNJG(3053) »sNBNJG(30) »BNL(3) s8NLG(3)sBNS(3)> 
e  RMEX(838)5UVEC(8) »CVEC (8) »ARR( 99) .DLNJ(30,3); 

e ERRPSI(1500) 


STDMU 
CHEMU 
TMP 
SPEiGs 
VNU 
QK Z 
BNJ 
BNJG 
NBNJG 
BNL 
BNS 
BNLG 
RME X 
UVEGSGVEC 
ARR 
DLNJ 
BRR Pow 


DATA INPUT 


SiO-StATE CHEMICAL POTS OF S29ECTIES 
CURRENTS CHEMICAL POTS OF SPECIES 

TEMPS “OF "AVAITEABLE STO=STATE CHEMICAL POTS 
SPECIES *LIst 

STOLCH FCOEFRICTENTS 

TORAL §# «OF *COMPONENT MOLECULES “PRESENT 
#20F FMOCECUPES SOP PEACH SPECIES 

ENITEAL GUESSES AT MOLE NGOS 

Ser VAGrs OF SPECTES 

CURRENT TFOTAL MOLE NOS 

SUM _ OP CURRENT SPECIES MOLE NGS 
GUESSEO TSTAE MOLE NGS 

JACOBIAN MATRIX 

VECTORS 

INPUT ARRAY=FOR IBM SSP ROUTINE "“GELG* 
ENCREMENTAL “E0GS, GF “SPECIES MOLE NOS 
ERROR FUNCTION ARRAY 


READ 251.5 MSPEC. LCGOMPsNPHsNGESs NTEUPsNTMPsKDOT 
READ 2523 “TTAs *I1TZs NITER 
READ 252s (SPECS(I1)+I=1+5,MSPEC) 


READ 253 


y COVNUC Ts J) s92=P, COMP). t=1.MSPEC) 


READ 2554((NBNJG(N) >» ( BNIG(NsL) sL=1+NPH) ) »N=1s9NGES) 
READ 270s (BNLG(LFA)»s LFA=1 +NPH) 

READ 2565 (QKZ(1)+sI=1,LCOMP) 

READ 2575 (TEMP(I).I=1>+NTEMP) 

READ 2575 (TMP(1)+sI=1+NTMP) 

READ 258s R»PRES 

READ 2595 TOL »ERFLIMsXJLIM 

READ 2615 (STDMU(CIsJ) sJ=1sNTEMP) sI=15MSPEC) 


MSEC 
Een M > 
NPH 
NGES 
NTEMP 
NTMP 
KDT 
ITAsiITZ 
NITER 
Wet 
ERE aM 
XJLIM 
R 

PRES 


OFT SPECTES 

OF COMPONENTS 

OF PHASES (MAX=3) 

OF GUESSED SPECIES MOLE NUS 

OF TEMPS OF INTERES? 

OFRMTEMPS GF AVAILABLE STP=STATE CHEM POTS 
CODE, FOR “IDENTITY GR CIBERWISS Gr ATeEMr4 ss V1 Me? 
RLiRS® SAND CAST TEMPS UF INTEREST 
MAXPNGVGOFPITERATTONSSSPECIFIED (< 1500) 
ERROR TOLERANCE FCR *GELG* 

MINIMUM VALUE FOR ERROR PARAMETER 

MINIMUM ALLOWABLE MOLE FRACTION 

UNIVERSAL GAS CONSTANT 

PRESSURE 


a a a: a 


*% 


DO 2 N=1sNGES 
HQ) 3! LFA=1>»NPH 


ads ‘Se < gateme 
204 2y0M JATOT THSRAUD 
20% 3J0M 2429582 THESSUo 30 ue 


*aJ80" SKTTLOR Sze Mat, ABA yAarA TUSHT 
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(2a Seu, f=0e(OMGO Je Tete IL eT UAV) ) « 


(2RDU fas LUHOM s t= 20 eMPOLMG 1 c (MH )OLMBH)) 22S 


(H4r~ f=as «(ASD OIME? 
(M09 Je t=Tel TISAS) «© 


eOTS 


(GMS THs l= 601 )9MS3T) 2 TSS 
(QM Met=ete(LI9RNT) «TSS 
ZaAG.P «62S 


MiJLA eMItS3FSe JOT «O&S° 


[9B92Net=1e (SM STH ee TSaLe (wet PUMOTS) » 


@ainsse 3 
2ThSHMO9MO2 FO & 
(2=KAM) 2324H9 40 © 
204 43JOM @2139=92 aseeauD 370 § 
T2geSTHMi 3O e9m3T 30 & 

ETCN END STATS-GTS SJBAUIAVA 30 2968T 30 
‘ay7? ,*4M9T* FO =2iWwssHTo AO VFITHIG) O39, 3005 
fests7Hl AD Samat Teau GHA T2n14 

{0922 >) GPLaresee BuGtTaesyI AS OW KAM 
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WOTTIAR= S40M 2IGAWDIJA MUMIVIM 

TWAT2HOD 245 JASRAVINU 

BaAvVzesAaG 


2204sr=K Ss 0a 


23926 32 


HOH, IeA5 of GO” 


veuuuy 
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vey 
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vuu 


ai 


vuuuvue 


wuuyY 


ere 


NN=NBNJG(N) 
BNJ (NN>sLFA) 
CONT INUE 
CONTINUE 
DO SepeSA=1sNPH 2 
3 BNL(LFA)=BNLG(LFA) 
DO 5 LFA =15NPH 
BNS(LFA) = 0 
DO 4 J =1sMSPEC 
4 BNS(LFA) 
S CONTINUE 
LNP=LCOMP+NPH 
PRLOG = ALOG(PRES) 
De (240 etTeENe =-1TAs,1 bz 
GORGW G10 +6 isnt 


BNJG(N»sLFA) 


No 


BNS(LFA) + BNJ(JsLFA) 


6 ehRCTEMPCITEMP) .GE. TMP(1) . AND. 


e TEMPC(CITEMP) «LE. TMPCNTMP) ) 
PRINT 322s ITEMPsTEMPCITEMP) 
GOTO" 250 

% DOs IB=1,.5NTMP 
IF (TMPCIB) «GT. TEMPCITEMP) ) 
8 CONTINUE 


9 FR = (TEMPCITEMP)-TMPCIB—-1 ))/(TMP( 13) -TMPCIB—1 ) ) 


1G CONTINUE 
IHOLD = 0O 
DEFT SSeiI FER = 1sNITER 


COMPUTE THE CHEMICAL POTENTIALS 
14 CONTINUE 
DO 35 LFA=1,NPH 
DO 30 J=1.MSPEC 
XJ =e ABSCENIC JsLFA) / BNS(LFA)) 
XKLOGy= ALOG(XJLIM) 


GOT93 


GOTOG 


lin ss) ms OG MleLM) ~ XLOG = ~ALOGCXI) 


GOTO (25,20) »KOT 


20 USTD = FRXSTDMU(J,1IB)+(1-FR)*STOMU(J,I1B-1) 


GOTO 29 
5 USTD = STDOMUC(JsITEMP) 
Q9 CHEMU( J) = USTD + PRLOG + XLOG 


0 CONTINUE 
5 CONTINUE 


TEST FCR CGNVERGENCE 
(A) EVALUATE ERROR PARAMETER 
40 SUMQ=0. 

ERF = 3*MSPEC*ERFLIM 

DO 43 K=1,LCOMP 

SUMB=C0. 

DO 42 LFA=1sNPH 

SUMA=0. 

DO 41 J=1.MSPEC 
41 SUMA=SUMA4VNU(J5K) *BNJ(JsLFA) 
42 SUMB=SUMB+SUMA 


43 SUMQ=SUMQ + (QKZ(K) - SUMB) **2 
SUMN=0. 
DO 45 LFA=1»NPH 

45 SUMN=SUMN+(BNLCLFA) - BNSC(LFA)) 
SUMG=0. 


DO 47 J=1sMSPEC 
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DO 46 K=1.LCOMpP 

SUMB=SUMB + VNU(JsK) *CHEMU (K) 

SUMG = SUMG + (CHEMU(J) - SUMB) **2 
ERRPSICITER) = .5* (SUMQ + SUMN + SUMG) 
Te’ CERRPSTCITERJ4LGE.AERF) GOTO 54 


EVALUATE MOLE FRACTION CORRECTIONS 
TeeGiterP EG. tT) GAGOKTO S4 
DO 49 LFA =1,.NPH 
SCOUNT =F 0 
DO 48 J =1,MSPEC 
eS =S6N0C9, @EAIL 7 SNS(UFA) 
DXUN=SxKUN+* DENI TISLEA) 
TE CABSCOXS) GE. 0.56-05}) GO TG 48 
JCGOUNT = JCOUNT + 1 
CONTINUE 
Te CICOUNT site MSPEC) GO To 54 
CONTINUE 
PRINT 30C,. TEMP(ITEMP) 
PRINT 259%SUMQsSUMN.».SUMG 
PRINT 311 
OO SSSNR*®=14,8NP 
PRINT 3215 (RMEX(NR»NC) SNC =15LNP),s CVECC(NR) 
PRINT 311 
PRINT 271+ (CHEMU(J),/ =1eMSPEC) 
PRINT 310 
PRINT 2709 (BNS(LFA)»sLFA =15NPH) 
GO TO 196 


C- SETTUPRNMATRIX *RMEX* 


54 


89 


85 


90 
Cc 


OO7T70 K=1.LCOMP 

DO 65 I=1.,LCOMP 

BRT =O, 

DO 60 LFA=1,NPH 

SUM=0. 

06°55 =J=1.MSPEC 
ADD=VNU(Js1) *VNUCJsK) * BNI (JsLFA) 
SUM=SUM+ADD 

CONTINUE 

BKT=BKT + SUM 
RMEX(Ks1)=BKT 

CONT INUE 

DG 90 LFA=1.+NPH 

LCF = LCOMP+4+LFA 

D9 85 K=1+LCOMP 

SUM=0. 

DO 80 J=1.MSPEC 
ADD=VNU(J 9K) * BNJ(JsLFA) 


SUM=SUM + ADD 

RMEX(KsLCF) = SUM 

RMEX(LCFsK) = SUM 

CONTINUE 

RPMEX(LCEsLCF) = BNS(LFA) - BNL(LFA) 
CONTINUE 
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DO 140 K=1,LCOMP 
SUM=0. 
DO 135 LFA=1.NPH 
SUMA=0.6 
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130 
135 


1490 


145 


Esc 


SUMB=0. 

DO 130 J=1.MSPEC 
ADD=VNU (CJ eK) *BNJ(J,LFA) 
SUM A=SUMA#FADD 
SUMB=SUMB+CHEMU( J) * ADD 


CONTINUE 

SUM=SUM + SUMB =- SUMA 
CVEC(K) = SUM + QKZ(K) 
CONT INUE 

DO 150 LFA=1.NPH 
GBFN=0. 


DO 145 J=1,MSPEC 

GBFN=GBFN + BNJ(JsLFA) *CHEMU(J) 
CONTINUE 

CVEC(LCOMP4+LFA) = GBFN -—- BNS(LFA) 
CONTINUE 


SULVE (FOR UVEC 
REARRANGE DATA FOR "*GELG® 
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161 


DO LeorisaieL Ne 
JJ=(J-1) * LNP 
DG MWSS WW=1.LNP 


ARR(JJ+1I) = RMEX( 19d) 
UVEGCH)+t= CVECTL) 
CONT INUE 

CONTINUE 

PERN= 0 


CALL GELG (UVECsARRsLNPsl1,TOLsIER) 
te CER MSEG 590.) GOTO 161 

PRINT 2ols (LER 

GUT0') 250 


LCF=LCOMP+LFA 
DO 180 LFA=15,NPH 
FACT = 1. 
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162 


166 


170 


UI=UVEC(LCF) 

UAB=ABS(UTI) 

TR(GAS WE. 6200) GO TO 162 

PACT 92.10 74 LUA 

UZ TS=it2 258 UI 

DO 175 J=1.MSPEC 

SUMA=0 6 

Mt ENJUIsLEA) 7 BNSCLFEA) 
TECADS(KI) .GE. NUL IM). GO 10466 
KOS OCI IM 


BRLEGUISUIEA i= xXJUDM * BNSC(LFA) 
XJL=ALOG( XJ) 
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SOerre (1=7 <L.CcOMP 
SUMA=SUMA4VNU( Jo 1) *UVECCTI) 

te) Oe a UI + SUMA - CHEMU(CJ) 
DENG TERA a= IL J 

TRG SFG a2.) «~GO FOL 17s 
remepits efie.2Uz) GG TO 175 

RATIO T= (ABS CUZ +7 40S) 

TEIRCRATIC (SLT SEACH) IEAGY SE SR ATI0 
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LS 
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NOG 


178 
BTS 


180 


195 


CONTINUE 


APPLY CCMPUTED CORRECTIONS 

BNS CLEFA) 70's 

DO 11.79 *J=1,.5MSPEC 

ADD = FACT * DLNJ(JsLFA) * BNU(JIsLFA) 
BNIIGCJsLRA bh és BNJ(J,LFA) + ADD 
PERUBNIV Is UPA ), "ei, 107s)! . TEND CISUERA Jace 
BNS(LFA) = BNSC(LFA) + BNJ(JsLFA) 
CONTINUE 

ADD = FACT * UI * BNL(LFA) 

BNLCLFA) = BNL(LFA) + ADD 

IF *CBNIECISEA ) Asi. OO. BNL CLFA) = BNS(LFA) 
CONTINUE 

PRINT 311 

CONTINUE 

GO TO 50 


OUTPUT 


196 


133 
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230 
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259 
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CONT INUE 
PRINT 320 
PRINT 311 
DO 198 I=1,ITERs5S 
IA = I-11 
1B IA + 4 
LG = tet. A 
PRINT 301.5 IAs(CERRPSI(J) oJ) sJ=1513)sERRPSICIC) 
PRINT 312 
IF (ITEMP .EQ. ITA) PRINT 309 
PRINT 3135(L sL=1+5NPH) 
PRINT 311 
PRINT 302s((SPECS(J) + (BNJ(JeLFA) sLFA=1+NPH))> 
J=1,.MSPEC) 


PRINT 310 

PRINT 3033 (BNL(LFA) »LFA=1 sNPH) 
PRINT 324 

PRINT 311 

DO 230 L=1+NPH 

PRINT 323 5L 

PRINT 310 

DO 225 J=1+MSPEC 

RJ = BNICJ st.) 7EBNECL) 
PRINT 302s SPECS(J)sXJL 
CONTINUE 

CONTINUE 

STOP 


FORMAT(1015) 

FORMAT(10(1X3A4) ) 

FORMAT(C4F522) 

FORMAT (8(15sF5e1)) 

FORMAT (4F19.3) 

FORMAT(SF10.1) 

FORMAT(2F10.3) 

FORMAT(8E10.2) 

FORMAT (5F 1563) 

FORMAT(3F10.3) 

FORMAT(5F10.3) 

CORMAT (91 %s 1Xs/7771SC 8," TEMP = 8st OuleeXs * DEG RR?) 
SORMALO 110.41E1 Seca lO} cl os3) 
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FORMAT( f'=-',1X) 
FORMATC® %51X) 
FORMAT (*1%,1X) 
FORMATC///T41 s*SPECIES'.TS7s*MOLE NOS*s 
J7#T4&95 "PHASES 11052120) 

FORMATC//TS0s"* ERROR PARAMETER TRACE!®) 
FORMAL VC TTAMs GP1L0w Ss 10.3) 
FORMAT CTLs 1X77 /7 554" TENP C4126 *) =) *sFSe ls SXs 

"OUTSIDE RANGE OF TMP) 
FORMAT. -C'=" sTSUs\* PHASE! . 15) 
FORMAT (///T44s"*MOLE FRACTIONS!) 
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APPENDIX C 
GRAY GAS APPROXIMATIONS TO REAL GASES 
C.1 Termination of the Gray Gas Computations 


The computation of a mixed grayt+l-clear-gas approximation 
to a real gas is, in essence, approached through the pro- 
Gréssive “extraction of the Contributions of successive 
gray gases to the radiative properties of the real gas. 
This is effected by performing the computations within the 
beam length ranges in which the different gray gases ef- 
fectively control all variation of the radiative properties 


of the overall gas mass. 


As suggested by the discussion in Chapter IV, the order of 
'extraction' of the gray gases is that of increasing ab- 
SOrption Coctficients, K.- Successive Ke are at least an 
order of magnitude different from one another, and so are 
the one : 

The computation of the component gray gases is to be ter- 
minated if the estimated maximum contribution of any as 
yet uncomputed gray gases (i.e., its weighting factor) is 
judged “anstgnzticant. Thus “the specification of this dis-— 


criminatory level of significant contribution annihilates 
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any need for an a priori statement of m, the total number 
of gray gases required. 

the extraction of each gray gas is carried out within the 
range of beam lengths in which the variation of its con- 
tribution determines the variation of the overall gas 
property. The upper limit of this range of beam lengths 
can be estimated for each gray gas except the very first 
to be ‘computed. At this upper limit, the contribution of 
the gray gas is qust beginning sto be Ansigniticant,- for 


example, at 5 per cent. of itsS maximum value; 


i. Gey, fae e los eee (om) 
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ont F max, i i 

But Rive Oyhsen 
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ererore max,i Clee 


As a check on the decision to terminate the computation 
of the gray gases, it may be worthwhile determining the 

: : hich the contribution of the un- 
beam length al at whi 
computed (ith) gray gas is at the minimum allowable level 


(0.0005, say) 
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Thus the decision to terminate would be vindicated if 
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C.2 The Listing of GRAYGAS 


The graygas approximation program, GRAYGAS, employing the 
technique set out in Chapter IV and this Appendix, is 
ered anetnis) section. lb 2S a, FORTRAN 1V program that 
was successfully run on an IBM System 360/67 computer 


instalied at the: Computing Services of the University of 


Alberta. 


147 


ar) 
eS 
SAQVARD to padteit ef? 4.9 
| 7 
ot pniyolqms ,2A0YAZD ,wtexpoxq noltsmixncaggs espystp ae 
et \xibnsagh aif? Bas VI SedqsH9 ab duo tee supinioad - 
‘Ss 


jad mexporq VI “ARIMOT s at ¢I .nofsose efdd mb Bedale 


secugmon Ta\0ae moxey2 MET as no nx livteessoue esw Oe 


2 ae 


26 yiteisvinw! ot to seotvise pricuqmol edd gs beliasent _—_ 
7 


@ ta ©) OY @ OO AO) © 


(O(a) 


QLOiG) OF OL OLOZOLOLOCOLOLCOLOLOLOFOLOLOLOFOLO: 


118 


GRAYGAS 

-COMPUTES A MULTIPLE-GRPAY-9SLUS=-CLEAR-GAS APPROXIMATION 
TO BURNED GASES IN A COMBUSTION CHAMBER OF GIVEN 
GEOM USING AS MANY GRAY GASES AS REQD 
(MAXIMUM 5) 

-FITS A TEMPERATUSE POLYNOMIAL TO EACH GRAY GAS WT-ING 
FACTOR OVER A GIVEN TEMPERATURE RANGE 


DIMENSIGN ,EC(30,20),EW( 30.20) s26(5530) ,AGI (5530); 


e GAM(5Ss5} «GMA(5e5) ,DELTA( 30520) sDELTB( 30.20). 
ss PLC(30).PLW(39)-.PWR(10),PCWL(15)-> 

° TGP (20) +TG(30)sTDCS); 

“ BML(20)+SIGMA( 20) ,SLK(5) «EGAS( 20) +VGI(20})>, 

. AGZ(€390) sAGR( 200) »GAR(39)s 

< NPVT CS) sELUCS) Ree (S)s. VECCS) .VECSUS) 


COMMON /XFE/1IB»+JBsEAsEB 


2 /POL/NTGsNGRAY »MP 
ECsEW EMISSIVITY (EM) TABLES FCR CG2sH20 
EG GAS EM VS BEAM LENGTH (B-L) VS TEMPERATURE (TEMP) 
AGI »AGZ EM WT-ING FACTORS 
GAMs,GMA MATRICES 
DELTAsDELTB EM DATA CORRECTION TABLES 
PLC>sPLW PL-VALUES FOR CO2sH20 EM*S 
PWR PRESS-RATIOS FOR AVAILABLE H2Q 
PRESSURE-BROADENING CORRECT IONS 
PCWL PL-VALUES FOR EM CORRECTIGNS 
TGP TEMPS FOR AVAILABLE EM DATA 
TG GAS TEMPS OF INTEREST 
TO TEMPS OF AVAILABLE EM CORRECTION DATA 
BML BEAM LENGTHS (B-L) 
SIGMA B-L/MEAN B-L RATIOS 
SLK GRAY GAS EXTNCTN COEFS 
EGAS GAS EM VS B-L TABLE AT MID-RANGE TEMP 
VGI VARIABLES 
AGRsGARsNPVT 9 
VEC.VECS VARIABLES 
FLU,RCF WORKSPACE (FOR CSLNIL) 
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C DATA INPUT 


c 


(a) (@) CY fay (@) (es ial) Tal (al) OO) (ek falta) CY (a) (Cr @ faba 


ox} 1) (Oj) Cem {hh Ce 


NTG 
NGR 


Se-=" 1 


FORMAT (515) 
FORMAT (F10.35,C9FS5.3) 


FORMAT 


(SE1S.6) 


FORMAT (2F10.3) 
FORMAT (SF10.3) 
FORMAT (6F10.5) 


READ(S5s1) NTGsNGRsNSSGoNQ 


READ(S.1) 
READ(551) 
READ(536) 
READ(5,56) 
READ(554) 
READ(S,4) 
READ(554) 
READ (5,2) 
READ( 536) 
READ(5.6) 
READ(555) 
READ(5,5) 
READ(5,6) 
READ(Ss2) 
READ(536) 
READ(536) 
READ(5,6) 
READ(5,57) 
READ(S557) 


NSSG 
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NTGP 
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FOr GAS TEMPS OF INTEREST 


MAX # OF GRAY GASES 
# OF B=-L/MEAN B-LO RATIOS CINSTDERED 


MUN 2 Or TERMS IN =o 1R-LINE—— 1a) 


POOR TEMPS FOR GIVEN COZ 
# OF PL-VALUES FOR CO2sH20 EM DATA 


MAX DEGREE OF TEMP POLYNOMIAL 
APPROX TO THE WT-ING FACTORS 
ENTRY # OF MINsMAX TG OF INTEREST 


DIAMsAXISs 
DND »XND 
PDBsPERM, 
PLI 


H2O EM*S 


DIMENSIONS OF COMB CHAMBER 


MIN FRACNL AGREEMENT BETWEEN: 


CCFACs CWFAC 


PHI 


=—PREDICTED €& ACTUAL VALUES AT MAX B=L. OF 
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FOS Sui 
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COZ. 20 EMSs 


INVERSE OF EQUIVALENCE RATIO 


GET MIDDLE RANGE TEMP 
RANGE =*ONFA4t NE) 7 2 
MIDR = RANGE 

TGMR = TG(MIDR) 


FR = 


TPAGFR SLT 2 EO? PPR 
DEN = 1 + 23.8*PHI 
ZEA 0.0 


(TG-MIDRANGE) 
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PRC = 3xFZ 
PRW = 4XFZ 


PTOT = PRC + PRW 


PR = PRW/PTOT 


PRE Oo 
GALL FIND (PWRs6sPRsPRFs001) 
JBS = IB 


GET JHE? MEAN B-L 


XLM = BEAMCAXISsDIAMsXNDsOND) 


PRINT 3259 XLM 


PRINT 360s TGMR 
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TFR = O. 
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PRINT 29660) TGMR 


GO TO 200 
1Su.JBX. = 1B 


DO 45 ISIG=1s,NSSG 


JB = JBX 


BM = SIGMACISIG) * XLM 


BML(ISIG) = 8M 


PC = PRC * BM 
PW = PRW * BM 
Figk= §PGEtePw 


EGX = OO. 
KD = 1 
KE = 1 
PFR = 0. 
KODE = 0 


CALL FIND (PLCsNPC, PC»PFRsKODEs1) 


lie CKODE VNB. 
PRENT 295 


xP. = pc 
UPA PLCC) 
Gorm? 27 
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GO TO 20 
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FAC = CCFAC 


GOTU 8235 
25 Ke = 82 

PFR = O« 

KODE = O 


CALL FIND (PLWsNPWsPWsPFRsKODEs1) 


IF (KODE .NE. 

PRINT 300 

XP = PW 

UP = PLW(1) 
27 PRINT 276% XP 

PFR = XP/UP 

KE 2 

GOTO (20.305 
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30 CALL EINT (EWsTFR»sKE) 


FAC = CWFAC 
35 ADD = Oe 
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EGX = EGX + AOD 
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KODE = 0 : 
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CALL EINT(CDELTA»sPRF,.01) 
DELA = PLF*EB+(1-PLF) *EA 
GALt EINT(DELTB»PRF»01) 
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Tee GEGXeeiay DEL) DEL = 0.5*(EGX=-DEL) 
EGAS(ISIG) = EGX-DEL 

1Ge (RU asisT. BPeWkEG2).) ) edEb s=<0% 
CONTINUE 

PRINT 320 

PRINT 2203(8ML(1)+,I=1sNSSG) 

PRINT 220% (EGAS( 1), I=1+sNSSG) 


CAL CUBRAT ESRHEIEXTNEIN COEFS 


SO 


Sa 


NSr= i1 

MUg=y1 

NU = 2 

AGSUM = EGAS(NSSG/3) 


AGSTP = 0.20*(EGAS(NSSG) -AGSUM) 
PER = PERM 

PEL = ALOG1O(PER) 

PDL = ALOG10( PDB) 

PLLN = ALOG(PLI) 

AGOs4= 74-0 

NGRAY = NGR 

SMX 9=nN0—9 

DO 70 IGRAY =1,+.NGRAY 

ie GCHGRAY _siGie 1) GOTO 51 

a = 0 

KSe=40 

CALL FIND (EGASsNSSGsAGSUMsF sKSo2) 
PRINT 240 

PRINT 216s AGSUM 

PRINT 241 

Te GKSSeEQsS .OR. IB.LT.2) (GOTO 80 
LFe(SIGMACIB) .«LT. SMX) GOTO 80 
S$Oe=-810*s0e2 

SREF = 0.01*SMX/SD 

5Oess. L=1. NS 

SREF = SREF*SD 

Sp=_SREE 

fee(l 2eeQ<nh), S = 0 

BML(I) = S*XLM 


FR = 0. 

KS = 0 

EGX = 0. 

CALL FIND(SIGMAsNSSGsSsFRsKSo2) 


feat beecGtsrlivecx = EGAS(IB=1)*(EGAS(IB) 
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GAR(I) = EGX 

CONTINUE 

het (BMEGNS) }.GTLOBNEGe)) GOTO 54 
DO 58 1I=2.,NS 

J = NS-I+2 . 

HB = BML(I) 

BML(I) = BML(J) 

BML( J) = HB 

HG = GAR(I) 

GAR(I) = GAR(CJ) 

GAR(J) = HG 

CONTINUE 

DOJ60 T=i1sins 

ARG = AGSUM - GAR(I) 

IF (IGRAY .EQG. 1) GOTO 56 
IGRA = IGRAY - 1 

DO SS JG=1.,IGRA 


ARG = ARG - ACGI(JGs1) *EXP(-SLK(JG) ¥BML(I)) 


CONTINUE 

TRE CARG) 52,52, 59 

te €iGRAY .~EG. 1) GOTO s0 
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SMX = BML(I-1)/XLM 


GUT St S52 
NGRAY = IGRA 
6GTG 71 


VGI(I) = ALOG10OCARG) 
AGR(I) = ARG 
CONTINUE 


SK=(VGI (NS-1)-VGI(NS-2) )/(BML(NS-1) -3ML(NS—2) ) 
VGQ = VGI(NS-1)+SK* (BML (NS) -BML(NS-1)) 


DB = VGQ-VGI(NS) 

IF (CIGRAY/NU) .GE. 1) GOTO 65 
IF (-ABS(0B)/PDL-1)65372572 
CONTINUE 

NSQ = NS-NQ+tl1 

00 6840=15NSQ 


A = 0. 
B= 0. 
C= Oe 
D= 0. 


DO 66 I=t>sNS 
BL = BML(T) 
VG = VGI(1) 
Aes A 42° BE 

B B + BL*BL 
Gh=nNG £aVG 
D 
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= D + BL*VG 
ONT INUE 


BOT = A*A-B¥ZL 

SLN = (A*C-D0*ZL)/BOT 

AGL = (A*D-8*C)/BOT 

ERRS= 0% 

HOté6é7 YVEESNS 

ERR = ERR+(AGL+SLN¥BML(J)-VGI(J)) *¥2 
ERR = SGRT(ERR/ZL) 

IF (-ERR/PEL-1)69968368 

CONT INUE 
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NOZ=ZE 

SLKCIGRAY) = ABS(SLN) ¥ALOG(10.0) — 
AGICIGRAY 51) = 10**AGL 

ira (E GRAYs EQ. 1) PRINT 340 
IFCIGRAY.EQ.1) PRINT 220. DBsERRsPDBsPER sPLI 
PRINT 306% IGRAY 

PRINT 220s (BML(1I),I=1.NS) 

PRINT 220+ (AGR(I)s I=1.NS) 

PRINT 308% tL» BML(L) 

PRINT 307sSLK(IGRAY) »sAGICIGRAYs1) 
AGO = AGO-AGI(CIGRAY 51) 

XLI = -PLLN/SLKCIGRAY) 

XLJ = EMIN/AGICIGRAYs1)/SLK(IGRAY) 
PP? (XU .LEs XLJ}) NGRAY = IGRAY 
SMX = XLI/XLM 

PRINT 230% XLI 

CONTINUE 

CONT INUE 

PRINT 2255 NGRAY 

NGC = 0 

GO TO 85 

GOTO (79s78+76)s MU 

AGSUM = AGSUM-DB*AGSTP/(DB-DA) 
GOTO" “7:7 

TRL UCD A=DZ) =( CE=DA)MSGT. 0) GOTO 78 
XA = 25% (DB-DZ)/(-DZ+2*DA-DB) 

XB (DB-DA) *(CA-D0Z) 

XC = (AGSTP+XB/AGSTP) *XA 

AGSUM = AGSUM-AGSTP4XC 

Nw@o 1 

GOTO 81 

TEA CDAZDE SET. ©} GOTO 73 

pz = DA 

DAY =" DS 

IF (MU «LT. 3) MU = MUF1 

AGSUM = AGSUM + AGSTP 

CONTINUE 

TRINMLAGY. EG. 50) GOTO 82 

PRINT 335,% DB 

iF “{UABC LEG. 562) GOTC 82 

PRINT 330 

PRINT 2155 (BML(1I)sTI=15NS) 

PRINT 2205 (AGR( I), I=1»NS) 

Pe CAGSUM .LE. 1.) GO TO S90 

GO TO 200 
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PRINT 240 

PRINT 310 

PRINT 311 

ES = 0% 

O06 SOPM=14NSSG 

BML(I) = SIGMAC(I)*XLM 
EGC = 0% 

DO 89 J=1sNGRAY 

PJI = SLK(J)*BML(T) 
AD AGI(Js1) 

SUR— = Ol. 

feu (Pulte LE. 175.0) SU = ADXEXP(-?2 JI) 
EGC = EGC+AD-SU 
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CONT INUE 

PRINT 315+sBMUC1I).EGAS(1I),EGC 
E = 100*(EGC/EGAS(1)-1) 
PRINT 312.5 E 

CONTINUE : 
PRINT 245 

TH ONGG .EQ. 1) GOTO 134 
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Eat 


BOL RSG PTG=PpeNntcG 

FR = (TG(ITG)-TO(1))7(T0(2)-TDO1)) 
BECGFR sb re#O)-4 sER=0'S 

TFR = O 

KODES= 0 

CALL FIND (TGPsNTGP,TGC(ITG) sTFRsKODEs1) 
if CKODEGSNE. ESN-8GG TO 94 
PRINT 305 

PRINT 276% TG(ITG) 

GO TO 200 

JBX = IB 

BMX = SIGMA(NSSG) *XLM 

DO 125 1I=1»,NGRAY 

BMN = —ALOG(PLI)/SLK(1) 

BML(I) = 0.5*(3”X+BMN) 

BMX = BMN 

IF (NGRAY .EQ. 1) BML(I) = XLM 
JO f= 78x 

PCL = PRC*BML(I) 

PwWL = PRW*BML(1) 

PL = PWL+PCL 


EGX = Oe. 
KD = 1 

PERT TOs 
KODE = 0 


CALL FIND (PLCsNPCsPCL»>PFRsKODEs1) 
TFi CKGDE=.NES59e) GO TO 100 
PRINT 295 


tes PCL 
UP =sPi2C C1} 
GOMTO 4125 


GAULETEUNT <CEGsTAR.KE) 
FAC = CCFAC 


Gowen saial 
KD f= y2 

PFR = O.« 
KODE = 0 


CALL FIND (PLWsNPWs PWLsPFRsKODEs1) 
BEiGKODE).NE. 9) GO TO 120 
PRINT 300 

<PHSRP WE 

UP = PLW(1) 

PRINT 276% XP 

PFR = XP/UP 

KE = 2 

GOera Grucds 20) sicko 

CALL EINT (EWsTFRsKE) 

FAG = .CWFAC 
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122 2GX = EGX + AOD 125 
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PLF = 0 
KODET=30 
CALE FINDO(PCWLo1C »sPL.PLFeKODE o1) 
IRPCKODETANE..€99)8 GOUTOTHI2s 
tsp =S10 
PUFT=F1 20 
les Yst=i19es 
124 CALL EINT(DELTAsPRF,01 ) 
DELA = PLF*=B+(1—-PLF)x*EA 
CALL EINT(DELTBsPRF,01) 


DELB = PLFX*ER+(1-PLF)*EA 
DEL = DELB*FR+t(1-FR)*DELA 
LEDPCPEV SLT. VeCGWEC2)) DEL = 0. 
PPMCEGXME «CT. DEL) DEL = 0.5*(CEGX=DEL) 
EG (IsITG) = EGX-DEL 
125 CONTINUE 
130 CONTINUE 
HOwi Z1PT=1 ¢NTG 
PRINT 285s TG(I)> 
e (EGC Js), J=1sNGRAY) 
131 CONTINUE 
DOLISS I=l.NGRAY 
DOt1T3S2 °9=1,.NTGE 
TE IUEGOIs)) «CE. EGMIN) GOTO 132 
NGRAY = I-1 
PRINT 225% NGRAY 
NGGt=#1 
GOTO 85 
132 CONTINUE 
133 CONTINUE 
CALCULATE THE WT-ING FACTGRS 
FOR ALL TEMPS 
WITH CONST EXTNCTN COEF 
134 DO 140 I=1s,NGRAY 
DO 135 J=1sNGRAY 
PJI = SLK(J)*BMLC(T) 


SUL= 0. 

Te) (G00 cbs PTS. OV TSSU T= ExXPC=PJL) 
GAM(IsJ) = 1.0-SU 

GMA(I+J) = GAM(IsJ) 


135 CONTINUE 
140 CONTINUE 
PRINT 240 
PRINT 345 
DO 141 I=1.NGRAY 
147 PRINT 215+. (GAMCIsJ)> J=1sNGRAY) 
PRINT 241 
IF (NGRAY .EQGe. 1) GOTO 157 
TEST FOR ILL-CONDITIONING 
DO 144 I=1+s+NGRAY 
SUM = A 
DO 142 J=1+NGRAY 
142 SUM = SUM4+GMA(I3J)**2 
RMS = SQRT(SUM/NGRAY) 
DO 143 K=1sNGRAY 
143 GMA(IeK) = GMACIsK)/ZRMS 
144 CONTINUE 
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M = 0 126 
CALL CSLND (GMAsSsNGRAYsDTs Ms GMA) 
PRINT 309% DT s™ 


PRINT 250 
IF CABS(DT*10**M)-DTMIN) 146,146,155 


USE APPROPRIATE CSLIB SUBRTNE FOR ILL-COND' EQNS 


146 
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DO 152 J=1sNuEG 

DO 147 I=1,NGRAY 

VEC C1) =Secerss) 

TONG INEIGT< 1) “GOTO 149 


GAEG CSLNILCGAMs VECs VECSs NGRAY »FLUsRCFeNPVT) 
GOTO 1350 

CALL CSLNIS(CVECs VECSs NGRAY ) 

DO 151 K=1+sNGRAY 


AGG J }) = 7VEGS GK) 
CONTINUE 
GOTO 175 
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GOTO CIS7#s Sega. hrash7syisn NGRAY 


GMA(19s1) = 1./GAM(1 91) 

GOTO 159 

© Tir=40'5 

Mi S=i00 

CAMs CSLNI(GAMs5 sNGRAYsDT»MsGMA) 


PRINT 3095 OT»s™ 
PRINTS S9 

DO 162 J=1sNTG 
DG 161 T=1s,NGRAY 
AGI(IsJ) = O. 

DO 16C K=1s,NGRAY 


160 AGI(IsJdJ) = AGICIsJ) +GMACIsK) ¥EG(KsJ) 
161 CONTINUE 
162 CONTINUE 
SOTO, 7S 
USE APPROPRIATE SSPLIB SUBRTINE FOR WELL-COND EQNS 
L7S CAL XMOVE (GAMsGAR»sNGRAY »NGRAYs591) 
CALL XMOVE (EG sAGRsNGRAYsNTG930 91) 
PER w= TO 
GAG GELG(AGR+sGARsNGRAY »NTG30.5E-C6e IER) 


FEMGLER PSEGs,.0) GOTO 174 
PRINT 2065 IER 


GOTO 200 

17S CALL XMOVE (AGI sAGRsNGRAY sNTGos 3022) 
17S 400 "177 11 FG=15NLG 

AGZCITG) = 1.0 

DO 176 J=1+sNGRAY 
R76 AGZGULTG) = AGZCITG) «= AGLGJ» 1 TG) 
177 CONTINUE 
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PRINT 240 

PRINT 270 

PRINT 275s (SLKCIGRAY)sIGRAY =1,NGRAY ) 
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CONTIN 
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CONTIN 
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2554) 1GCTTG)+sAGZU TEN 
CAGICIGRAYsITG)»s IGRAY=. sNSRAY) 
VE 
POLY (TG sAGI » NTGs NGRAY » MP) 
UE 


( 110) 

(F10.1') 

(/10F13.5) 

(//TSSes*AGSUM = *,F10,5) 
(/10E13.5) 


(///T10s" TOTAL NOD OF GRAY GASES = *,12) 
(///T2s*BEAM LENGTH FOR MIN GRAY GAS °*, 

*EMEGONTRE=*. 67S. 5) 

Ct—a) 

(//) 

Cogs 
(T53s11H*¥*XRESULTS**///) 


(77T!STS (PHT) = SEG et) 

(/26Xs6HDIAM =sF5.2910Xs8HLENGTH =9F5.2) 
(*-", 5Xs*GRAY-GAS EXTNCTN COEFS®) 
(10X,5F20.6) 

(10X3s3F10.5) 

(////745Xs *GRAY-GAS EM WT-ING FACTORS!) 
(/T20s"GAS TEMP®) 

CUP (S0.t helt t=") 

(T229°*DEG- Rt) 

CEOS T414°O* 51 SX0*1% 613X082 1 3X0 4S*513SXs 
"4,15X_o"5!*) 

C/F27 o19T30+ 6614.5) 

(10X5***PLC TBL TOG NARROW**®) 
(10X.***PLW TBL TOO NARROW***®) 

(10Xs*** TG TOO LOW FOR TGP TBL **) 
(/7T2 s* GRAY GAS emer Be) | 
(//T2s*EXTNCTN COEF ="oF1C0.5s 

"EM WT—-ING FACTOR =',613.5) 

(//T2s"*STR LINE PLOT EFFECTIVELY STARTS ', 
ATCT #£ *s12% 

‘ CORR BEAM LENGTH = *#,£12.5) 

(7 * DETERMINANT "o9T25s"MATRIX 3 %5F10.25s 


(///T40s"*CHECK EG-VALUES's, 

7/7T2325s*BEAM LENGTH® »T50+"EG-OLD*" sT605 

*"EG-APPROX!?) 

(*4'.7T76s"ERRORIS/) 

(*4+",7T700F10.2) 

(TAS -EIZSSIER2SSs21SS5) 

(///T2s"A DISPLAY OF BEAM LENGTHS AND ', 
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APPENDIX D 


MEAN ABSORPTION COEFFICIENTS KR AND Kor AND FTFIELDX 


D.1 The Rosseland-mean absorption Coefficient, Ke 


The definition of Kn is given in different forms by Abu- 
Romia and Tien [1], Deissler [10] and Hottel and Sarofim 
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of the m wavebands, Ad; - If a mean value, K., is used for 


K in each waveband equation (D.1) becomes, 
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Which, by analogy with equations (D.1) and (D.2), readily 


leads to 
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This program is listed below. It is a FORTRAN IV program 
successtully run on an IBM System 360/67 anstalVation at 


the Computing Services of the University of Alberta. 
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GALL EQNS4(KSUB) 
PRINT 75 
PRINT 60+KSUB 
bEVCTIRESCTteel? GOTC 20 
FILL RHS VECTOR 
CALL RHSPS (KSUB) 
PRINT 80 
PRINT 60+KSUB 
SOLVE MATRIX EQN 
20 CALL SOLV6 CITRsERR»sKSUB) 
PRINT 85 
PRINT 50s KSUB 
PRINT 8&6 ITPR»ERR 
NGE=s (0 
DGOgM 22 411s NZ 
NEAR AEZNC AD). Gik: OF JE 46GCTOQ 22 
Doe CHUELNE «27 Oy GOT Dyee8 
Nie =) 1 
PRINT 125 
24 PRINT? 1305, a 
22 CONTINUE 
LESL ERR SLT we Ss OD) mGOTR; 25 
Like = ITR+1 
feaCidRyebi lke dd ERMJee GOTO 15 
PRINT 135 
PRINT 140 
PRINT 11435 (GI stTZNGld Jat =1 3 NZ) 
PRINT 98 
PRINT 145 
DOP 24, T=1%iNZ 
24 PRINT 1005 (EARCIsJ) »J=19NZ) 


SOTOnay 20.0 

CAL GULA TEMSURRACE FLUXES 

25 »CALL QFLX7 (KSUB) 
PRINT 90 


PRINT 605KSUB 
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PRINT 95 138 
PRINT 105 
DO 30 TI=13sNE 
PRINT 99,41 
30 PRINT 1005 CEAR( IsJ) » J—=15.NZ) 
PRINT 106 
DO 25 I=NE1.,NES 
PRINT 99,1 
35+ PRINT 100" CEARCIGU) JH1IGNZ) 
PRINT 1C7 
DO’ 40 I= NES1.NZ 
PRINT 99, I 
40 PRINT 1C0Os (EAR( I oJ) sJ=19NZ) 
PRINT 110 
PRINT 105 
PRINT 11245 (CIsTZNCI) )sI=1.NE) 
PRINT 106 
PRINT 1145 (CIeTZN( 1) )>TI=NE1.NES) 
PRINT 1C7 
PRINT 114+ (CCIsTZN(1) )s I=NES1 »NZ) 
PRINT 120 
PRINT 105 
PRINT 1153 (CisQFX( I) )sI=1,NE) 
PRINT 166 
PRINT 115+ (CIsQFX(1))sI=NE1.,NES) 
45 CALL EBLC8S (KSUB) 
PRINT 91 
PRINT 60,5 KSUB 
50 CONTINUE 
55 FORMAT (//T2s*ZONIR1I ?3*) 
60 FORMAT (*+*,T15s"SUBRTNE RC=" 512) 
65 FORMAT (//T2s"*NEIBR2e :*) 
TO BORMAT 777 2s'EXCHS :*) 
75 FORMAT (//T2s"EGNS4 3*) 
80 FORMAT (//T2s*RHSPS 3°*) 
85 FORMAT (//T2s"SOLV6 :*) 
86 FORMAT (//T1ICs*ITERATION #%5145T355"RMS ERROR =", 
oP RETSSS) 
GO FORMAT ™(//T2s*QELX7 :*) 
91 FORMAT (//T2,°EBLCS 3:*) 
95 FORMAT C€# 1" 5/47/7555 § ¥*¥*¥RESULTS#X*"S/T46 5 
- CLP ZEXCHANGE AREAS'/) 
98 FORMAT (*1%,1X) 
99 FORMAT (T2s12) 
100 FORMAT (/T2,10E£13.4) 
10S FORMAT(/T1,*"END ZONES") 
106 FORMAT(/TIs © WALL "ZONES® ) 
107 FORMAT(/T1s"GAS ZONES*) 
110 FORMAT eCl1'./7/7750s'(2) = ZONE TEMPERATURES *// 
7B T4S9°ZONE NO. 'oT70s "ZONE TEMP. "/T7Os1 IHC DEGREES R}) 
414 SORMAT (CT45.15.sT70sF9.2) 
115 FORMAT (T45515,1T60:E20.5) 
130) EORMATIVG! 1° .777/TSS0'(3) = ZONE FLUXES*// 
T455°ZONE NOs "sT705"ZONE FLUX DENSITY®/T70 sSH(BTUS >» 
9HSQ FT-HR)) 
125 FORMAT (///T2s* INVALID ZONE TEMPERATURES*/T2¢ 
~ PORT ZONE NOS ss * 
130 FORMAT “C/T2rT10) 
135 FORMAT (///T2s35HSPECIFIED NO OF ITERATIONS EXCEEDED 
~26He eo eEXECUTICNTERMINATING ) 
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140 FORMAT (/T2s3SHCURRENT TEMPS (DEG R) ARE AS FOLLOQW/ 


e 4T45s6HZONE #3T70s9HZONE TEMP) 


145 FORMAT (//T2s17FEXCH AREAS FOLLOW) 


200 STOP 
END 
SUBROUTINE ZONIR1 (KSUB) 
DIMENSION XMN(25)+RMN(25) 
CCMMON 7 AISRMNs» XMN 


e 7BISAXsDI SB2/DXsDR /B3B/NXGsNRG 
° /B4/NZsNESsNEsNSeNG /BS/NE1SNES1 
° 7C1ZRM 
GAL CULATE TOTAL «NO. GF ZONES 
NE = NRG+NRG 
NEl = NE+1 
NS = NXG 
NES = NE+NS 
NES1 = NES+1 
NG = NXG*NRG 
NZ NES+NG 


GENERATE «CGGRDS.~OF EACH ZONE 
(A) END ZONES 


RAD = DI/2 
DR = RAD/NRG 
NOTE: NRG IS TO BE SUCH THAT RM < DR 
PO LO; b=1 Nes 2 
Gig= gl 
RMN(I) = (X1/72-0.5)*DR 
10 RMN(I4+1) = RMNC(I) 
DO 15 I[=1.,.NEs2 


XMN(C I) = Of 
1S XMN(CI+1) = AX 
(B) WALL ZONES 
DO 62.0" T=NE1 s NES 
20 RMN(I) = RAD 
DX = AX/NXG 
A=) =DX 
B@ 25, 1= NE1>NES 
A = A+DX 
25 XMN(I) = A 
(Cc) GAS ZONES 


DO 35 I=1.NRG 
R = (1-1)*DR 
Agee DX 
K = NES+( 1-1) *NXG 
DO 30 J=1sNxXG 
A = A+tDX 
K = K+tl 
RMN(K) = R 
30 XMNC(K) = A 
35 CONTINUE 
EXHIBI TwtHE RESULTS OFPUAHIS, SUBRGUTIN: 
PRINT 40s NZsNEsNSsNG 
PRINT 45 
PRINT 50 
PRINT 555 (CIsXMNCI)sRMNCI)) »1T=15NE) 


PRINT 60 
PRINT S59 (C1 »XMNCI)»+RMNCT)) »IT=NE1,NES) 
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ao# (t.0<-SN2K) = C1)0MR. a” 
(to4#ms = (fede OF 
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«0 = {Lo“mMx 7 
“A = GietMMx Gr 
2av0s saw (6? 
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PRINT 65 140 

PRINT S50 (C1s+XMN(I),RMN(I)).I=NES1 NZ) 

KSUB = 0 

RETURN 
40 “FORMAW (///7T2s "TOTAL NO. OF ZONES =", 13s// 

e T1s"*BREAKDOWN FOLLOWS :*7 

oat 20s “END ZONES: 3 OST A/T 19, "WALL ZONES 2'6147 

a *T 206" GAS ZONES <$*,14) 
435 TPORMAT? 1C//T235"*ZONE NO.*sT62s "AXIAL COORD*®»T1025 

‘ "RADIAL COORD',s) 
50 FORMAT (/T1Cs*END 3$*) 
55 FORMAT (128,2F41.3) 
60 FORMAT (/T10+ "WALL :*) 
65 FORMAT (/T10s"GAS :*) 


END 

SUBROUTINE NEIBR2 (KSUB) 

DIMENSION RMN(25) + XMNC25) sNABGOR( 2565) 
CCMMON 7 AISRMN» XMN ZAS/NABOR 


e 7B2/DXsDR SB4/NZsNESSNE /BS/NEI2NES1 


INITIALISE NABOR 
DOt2J=1 45 
DOe AST=1'FNZ 
1 NABOR(IsJ) = 0 
2 CONTINUE 


(A) END ZONES 
DO 15' I=1,NE 
KO=? 1 
GOTOGNG 
Aels Ae2 NEIGHBOURING SURFACE ZONES 
3 DO 5° 9=13NES 
IF (RMN(J)-DR .NE.’ RMN(IT)) GOTO 5 
IF (XMN(J) eNEo XMN(I)CANDCXMN(J)4DX 2NE. XMNC(T)) GOTO 5 
Ke=0K+4 
NAGOR(I+K) = J 
S CONTINUE 
Ae3 NEIGHBOURING GAS ZONES 
6 DUNIONU-=NES1,NZ 
IF (RMN(CJ) eNEe’ RMNC(I)) GOTO 10 
IF (XMN(CJ) eNEe XMN( IT) eANDSXMN( JS) 4DX eNEe XMNCI)) 
eGOTO 10 
SIKT= KF1 
NABOR(I+K) = J 
10 CONTINUE 
NABOR(Is1) = K-1 
15 CONTINUE 
(B) WALL ZONES 
DD #35" 1=NE1s NES 
Kv 1 
GUT 26 
Bels Be2 NEIGHBOURING SURFACE ZONES 
16 DO 20° J=1-sNE 
IF (RMN(J)4+DR eNEe RMN(I)) GOTO ZO 
IE (XMNCJ) eNEe XMNCT)¢AND2XMNCJ)-DX ZNEe’ XMN(I)) GOTO 20 
K = K+1 
NABOR(IsK) = J 
20 CONTINUE 


DO 25 J=NE1+NES 
IF (XMN(J)-OX eNEe XMNCT) eANDSXMNCJ)4+DX .NES XMNCI)) GOTO 25 
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K = K+l 141 
NABOR(CIsK) = J 
25 CONTINUE 
B.3 NEIGHBOURING GAS ZONES 
26 BU? 30 "J=NES1,NZ 
IF (XMN(CJ) .NE. XMN(I)) GOTO 30 
IF (RMN(CJ)+0R .NE.’ RMN(I)) GOTO 30 
K = K+ti 
NABOR(I 5K) 
30 CONTINUE 
NABOR(I 91) 
35 CONTINUE 
(Cc) GAS ZONES 
DO SS I=NES1,NZ 
eta 


| 
= 


tt 
K 

{ 
~ 


Cork NEIGHBOURING END-ZONES 


DO 40 J=1.NE 
IF (XMN(J) «NE. XMNCI).ANDSXMN(J)-DX LNE!’ XMN(I)) GOTO 40 
IF (RMN(J) eNE.’ RMN(I)) GOTO 40 
K = Kt! 
NABOR(I.K) = J 
40 CONTINUE 


C €.2 NEIGHBGURING WALL ZONES 


DQ 45 J=NE1.NES 
IF (XMNCJ) eNE. XMN(I)) GOTO 45 
IF (RMN(J)-OR eNE.’ RMN(I)) GOTO 45 
K = K+1 
NABOR(I»K) = J 
45 CONTINUE 


CGC Cos SNEPGHEOURING GAS ZONES 


a) 


DG S0 HI=NESVseNZ 

IF (XMN(J) eEQ. XMN(I)) GOTO 46 

IF (XMN(J)-DX SNE’ XMNCI) -ANDSXMN( J) #DX SNE’ XMN(I)) G3TO 50 

IF (RMN(J) NE.’ RMN(I)) GOTO 506 

GOTG 49 
46 IF (RMN(J)-DR .NEe RMN(I)2ANDJZRMN(J)40R 2NE.’ RMNCI)) GOTO 50 
49 K = K+tl 

NABCR(I+K) 
50 CONTINUE 

NABOR(I51) 
55 CONTINUE 


{ 
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| 
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| 
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PRINT RESULTS OF THIS SUBROUTINE 


PRINT 70 
PRINT 75 
PRINT 80 \ 
Krav 
ci =e) 
PUV=FINE 
SOM Deveue b=1ls1U 
PRINT 85s IsNABCR(I51) 
JM = NASOR(I.1)+1 
IF (JM .EQ. 1) GOTO 60 
PRINT 905 (NABOR(I9J)sJ=2sJM) 
60 CONTINUE 
K = K+t1 
GOTO (61+62+63)>5 K 
61 PRINT 95 
Tie =eNeEd 


foo ((2)4MxX 934, 


TUNES 
GOTO S9 
62 §P- RI NieeLoo 
Il = NES1 
PU) = NZ 
GOTO S9 
63 KSUB = 90 
RETURN 


TO FORMATC*1*.7/7T405"ARRAY OF ZONES AND NEIGHB8DURS 


142 


FOLLOWS® ) 


75 FORMAT (77T142"ZONE*,T25s"* NO. GF*sTS0s*NEIGHBOURING ZONES#/ 


© 1149*NO.*sT25s "NEIGHBOURS! Ss) 
80 FORMAT (T1+s"END ZONES :*) 
85 FORMAT (1T1312eT30512) 
90 FORMAT ('+'%,7405;5110) 
95 FORMAT (Tls*WALL ZONES :!* ) 
100 FORMAT (T1+"*GAS ZONES ?:*) 
END 
SUBROUTINE EXCHS (KSUB.ITRsITN) 


CCMMON TAISRMNs XMN SAC/STZN ZAS/NABOR EAR s XDI 
° 7BISAXsDI SB2/DXsDR SBS/NXGsNRG SB4/NZsNES sNEsNSoNG 


e 7 Bo/NEVs NES 
e 7CISRM sHTCsES 


DIMENSION RMN(25)sXMN(25)sTZN(25) » 


e EAR(25525)sNABOR(2505)>5 
UAB (25. 5)5RSK(25; 5) > 
e XDI (646) »LW1(6) »LW2(6) 


INITIALIZE THE ARRAY "EARE 
Da 20° T=PaNZ 
DG 5) Js15NZ 
5 EAR (isd) 4—=20. 
10 CONTINUE 
VO = /AREA(DISO0I74) 
View— AREACDI+4DX) 
V2 1/ES-0.5 
V3 = AREA(2*ORsDX) 
HOL AREA( RM, RM) 
BOPSTeL= 7 1NZ 
Rie= ROSK(TZNCI) +1) 
Re =) ROSK(TZN(1) 22) 
RF = 0.75*R1 
M = NABGR(I,1)+1 
Se ethee.GT. 1) GOT 25 
15 PRINT 130% I 
20 KSUB = 1 
RETURN 
(1) GET A TWO-TERM REPRESENTATION OF 
RECIPROCAL GF UNIT—-AREA DIR-FLUX AREAS: 
25 If C1, GU. NES) GOTaeso 
TF aC hes.GTs NE) GOTO 40 
END ZONES 
360, B= 40x72 
DO 35 J=2>sM 
K = NABOR(IsJ) 
UAB(I+J-1) = V2 
RSK(I»sJ-1) = RF*S 
35 CONTINUE 
GOTO 157 
WALL ZONES 
40 Been DRAL 
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2=Mo5. 
S\5G = 6 OF 


IPO ONRG:) EQ... 1) B = OR 143 
DOR4S9JjJ=25M 
KM=eaNABBR CI ew) 


44 VABCIsJ-1) = -RF*B*¥B/DI+V24tRFEXRF/S(R2*RI*XDIS2) 
45 RSK(CI»sJ-1) = RF*B 

GOTO! 87 
GASY ZONES 


SOSD0655ed=25M 

K = NABOR(I oJ) 
(1.1) GAS-END 

Ba= Dx/2 

Cz= ve 

Ife (K. wISESY NE) GOTO 53 
(1.2) GAS-WALL 


B = DR/2 
IF (NRG .EQ. 1) B = DR 
C = -RFXB*¥B/DI4V24RFXRFE/(R2*¥RI*DIS2) 


De GK isle yg NES) GOTO 53 
(1.3) GAS-GAS 


B= OX 
Ce =r Ov 
IF (RMN(K) .EQ. RMN(I)) GOTOQ 53 
B = DOR 


S53. YARBCIEJ—)) 9c 
55 RSK(I+J-1) 
S57 CONTINUE 
(2) COMPUTE DIRECTED-FLUX AREAS 
(GRAY GAS CONTRIBUTIONS) 
DOE 7Of 151 5NZ 
AN = AREA(2*RMNC1I)4+DRsDR) 
M = NABOR(I.51) 
lif Cl GTS ONES)- GOTO) 6.1 
C2, Inisa G2 eae) GAS-SURF ACE 
Ae =. AN 
ie Giy wGT es sNEJA = Vi 
DO 60 J=15M 
K = NABOR(I,J5+1) 
60 EAR Clisk yee VAZAGUABC ta J) 4RSK Gls-1),) 
GOO’ 7 <7.0 
(2.3) GAS-GAS 
61 DO 65 J=1oM 
K = NABOR(I,J+1) 
Ait= FAN 
fe Ge. Guo besh GOTO ~62 
TE GiuseGhe NEA =. V1 
GOLO Mees 
62. tar= ak! 
ie Gs he O20 SAR 
be Cle ene Gayl) 1 (GOTO +65 
Av=7AXDX7DREVS 
ie GK esteb ppl) A = A-VS 
eo) BAR Gl aa) = (AZ CORSK Cle ) FUABCI«J))) 
790 CONTINUE 
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SESCUIR aor. 1) GOTO 121 

PRINT 131. ITR 

DOs7) f=keNZ 

M = NABOR( 161 )+1 

PRINT 132s Is(NABOR(IoJ)s J=23M) 


M= M-~l 
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PRINT 4533 
PRIN Te ttSs 
M = M+#l1 
PRINTS 3:3 
71 "CONT NUE 


GET DIRECT EXCH AREAS FOR SURFACE-SURFACE EXCH 


BY “ORS Ctra 


THENCE! TRANSFER MATRIX FOR TOT EXCH AREA CALCS 


DOGO 80 I=1 
lf <1 SEQ 
A2 = RMN( 
Al = A2 
DE Meters tore 
Ti = I+1 
DO 79 J=I 
B4 = RMN( 
B3 = B4 
Lee tele 
XJI = XMN 
XK = ABS( 
X13 = XK 
K = I 


Py CLWPNENE 16a SG SeSik Ss I 
seeCR'SKC 1s 5) 6 J=15) 


2 (EARCIsNABORCI oJ)) sJ=29M) 


CMPUTAT ION 


osNES 
e NES) GOTO 121 
1) 


e NE) Al = A1+DR 


1,NES 
J) 


- NE) B3 = B3+DR 
(J)-XMN(1T) 
XJI) 


Teen Test oO.) K = J 


PECK «GT 
X14 = XK 
Lee GxX Ia 
ti Cie. fe 
73 X23 = XK 
Pe CxXUI SS 
ita CIs 
74 X24 = XK 
K = J 
1G aed OF 3 9 ar 
BF GK «GF 


o NEY XisS X13-DX 


ie ere) EIR) Ae 
e NE eAND.J si Gilie NE ) X14 — 


GT's" os') GOTG 74 


CE. NEsAND.J «GT. NE) X23 


oto Oe) K = I 
e NE) X24 = X244+DX 


X14+DX 


EX = DIRECT(A1,B835X13)+DIRECT(A2 + B4 5X24) 
° -(DIRECT(A1sB4s X14) DIRECT (A29635X23)) 


XDI( IJ) 
TO MBDLCIs!) 
80 CONTINUE 


CORRECT sFOR 

rel .De Phaser 
IF (RMNCI 
F = HOL/A 
DO 124 J= 
Le Od! ¥2G 
LeU Re 
Pres GT 
AZ = RMN(¢ 
A2 = A3 
LECCE 
XIN = AX- 
X12 = XIJ 
LP OS St 
X13 = XIiJ 


= EX 
= EX 


END (EXIT) HOLE PRESENCE 
2sNEs2 
YAlGre "RMI GOTO 126 
REA(2*RMN( I)4+DR- OR) 

1 5NZ 

OY) eGoaTo @124 

Gholi)  GUTo i122 


* Nes? * coro * i122 
J) 

. NE) A2 = A2+DR 
XMN(CJ) 

J) NED « Kee Sale x 


A = DIRECT(RMs:A23X12)—-DIRECT(RM2A39X13) 


XOT Ole Jd 
XDICJs1) 


=) XDI( Is J}=A 
="XDI'CUs I) —A 
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hae 
184 TO. jean tance 
2 «836 FF 7 
tame 

AG+IA = tA im an 


agsea = 5S (34 484042 SE. 
(ToMMR=-CL AMX = TUX ' 
(TtRP2BA = AX 
ak EIx 
oa g 
t= ny yO s¥ss°NP-aD = 
MG-EUX = EX (aM .TOL AP SE 
aA = os? a a) 
Ev OTOD (oO .¥ie TURRET 
KOvetx = ore (3 .TO. 1eGHA, 9M «34. 1) AE 
MX = ESX ET = 
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9h = 88K OT 
et = 
I % a0 oTde TUR) VU 
Ate sx ASK LSA eTas BD FT . 
(42K, 6, SA)TOARIO+(CLAsE Ge 2A TOSAIO = xB ¥ 
((ES* Ea eSAPTIIAIOS (4 1% o6, LAD TOBAIO)=— Po 
aa = (t+1)TOX 
A3 = (ret )ITOK OF s ° 
SuuiTvod os . 


nou 


3J>MaeSzRe BUOH (TIX3S) Owe NOW TOSAROD 

SeIWMeS2t 28t OG ISt 

as! oTe> (wr ,To,. Clem) AI ‘ 
(AO.RO+( 1 IAM RSS TASAANJON = 8 

Sw.t=t «Sf oa 

ast Rahs 7 +02. 1 ad 7 
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GOT 
Al 
AJ 


Lie GuiR 


a) 


124 


EAR( I,J) xF 
EAR(Jo1)x*F 
BARI Cine) 
PAR (ust jie 
CONTINUE 
CONT INUE 


EAR(I»J)-AI 
EAR(JsI)-AJ. 


sil et 2} GOTO sel52 


DO 144 f=1,NES 


= 
Ir 
Di 
—— 
Up 
SU 


SU 


Vv 
GI 


(X 


1 
eGTe 


NE ) GO EG sl A2 


2*RMNCI)+D0R 
AREA(D1.4.DR) 


MN(T) 


eEQe. AXZAND.RMN(I ) 


-A*XES/(1-ES) 
DG 143 J=1.NES 
SU-XDI(I 9J) 
XDICIeI) = 


SU 


COMPUTE TOTAL EXCH AREAS 
(CLEAR GAS CONTRIBUTIONS) 
PRINT 136 
DO 14541=1.NES 


145 
ist 


146 
1S2 


147 


PRI 
CAL 


NT 
L 


134, 


PRINT 1385 
PRINT 1137 
DO 146 T=1,NES 
PRINT 134s; 
DOSESO I=1sNES 


pS 
= 
D1 
AOS 
iG 
AE 
AGI 
ite 
Jl 


Vv 
CI 


Ct 


(1 


1 
ofS o 


(XDICIsJ)sJ=1sNES) 


ole ts 


RM) A=A=HOL 


MINV(XOI eNESsDD>LW1 sLwe) 


DD 


(XDI (I sJ)sJ=1>sNES) 


NE) .GOTO-147 


2*RMNC1I)+4DR 
AREA(D1.DR) 


2 FQ. 


2) A = A-HOL 


AXES/(1-ES) 
AGO(TZNC1I)) 


2 EQ. 
I+1 


NES) GOTO,,1.50 


DO 149 J=J1+NES 


yal t= 
1G 
Le 
AF 
WE 


Vv 
(J 
(J 


1 
aILle6 
eEQe 


NE) A = AREA(2*RMN(J)+4+DRoDR) 


23) A = A-HOL 


AXES/(1-ES) 
-XDI(I,JS)*AFXAE 


EAR(IsJ) 


EAR(JoI1) 


149 CONTINUE 


150 


CK OX 
SELF EXCH AREAS: SURFACE ZONES 
END ZONES 
oO Ss AS lsNe 

AREA(2*RMN(1I)+4+0R+DR) 


Cc 
c 


EAR(I 41) 


je\ Bee 
Le 
ES& 


(X 


MN(I) 
ESxA 


AGI*WE 
AGO(TZNC J) ) *WE 


—-AGI*X(XDI(CI.1) *AE+ES) *AE 


eEQe AXSAND.LRMNCTI ) 


DO 154 J=15,NZ 


it 
ESF 


(J 


Ble 
ESF- 


c 


Ti asGOTD 154 
AR(IsJ) 


elvines 


RM) A=A-HOL 
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», Qe sy 
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Tal io 
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ive=a . 
T4f GTOSD (354 «TH— 97) BE 
AGeti veges = 1G J 
(20, 10)ARSa = A 
JOH-A = A {2 ,83, 1) a1 _*> 
(2a-1)\238A = SA TA! fs 
C{TIMRTIODA = IDA 
Cat oTOD {23M OB, 1) 31 
I4¢r 2 iu @ 
254+it=. C#f OO «<. 
IivVeza an 
{50 eH04(LIMMSOS AGRA = A (24 «.3d0 9 AL) a 
JOM-A 2 A 0S oOBa OF BP iy! 
{24-f)\eFea = FA 
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SOP4AGCL T)LGK= = By i 
34*17e 8 86=) = 6§(te bRAR 
Swei(LJHMSTIGOA = 41.4 3RAR 7 


alte 
3A+ (2343Ae(1. 1) 20xs"TOA- ~ tte tpAag oat 
235404 sibshiieeinaieh 4D: 


Cc 


154 CONTINUE 
TR ™CeESe sGe. 0.)  EARCIs1) = ESF 
155 CONTINUE 
WALL ZONES 
DO 160 I= NE1.NES 
Ae = VE 
ESF = A*xES 
DO 159) *J=Ts NZ 
Pee sea. i) GOTO 159 
ESF ="ESF-EAR(I. J) 
159 CONTINUE 
i Cesr  eGee Os) CARCI «1 ) = ESF 
160 CONTINUE 
SELF EXCH AREAS:GAS ZONES 
165 DO 170 I=NES1 NZ 
VK = 4*ROSK(TZN(1I).1)*OX*AREA(2*RMN(I)+OR30R) 
DO 1697 J=17,,NZ 
TP Coro. 1)? GOTO 169 
VK = VK-EAR(I4J) 
169 CONTINUE 
170 EARCTs1)" = VK 


Cx * 


@EG@) 


OD COR 


130 FORMAT (///T40o"** EXCH3 ERROR : ISOLATED ZONEs NOe'>s 
: I3,° ***) 

131 FORMAT (°1°,.T45s"EXCH3 : CHECK PRINTOUT®/T3>5 
. "ITERATION #*,13/7) 

132 FORMAT (/110.121,4113) 

13s PURMAT® (/T 253 SE13.4) 

1247 FORMAT  (772's51 0613.4} 

135 FORMAT (//T2s"*CLEAR GAS CONTRIBUTIONS :°*) 
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136 PORMAT = oss 125s" TRANSFER MATRIX FOR TOT EXCH AREA CALCS* 77) 


A "AREA CALCS*//) 
137 FORMAT (/S//T25s" INVERSE GF TRANSFER MATRIX#*/S/) 
128 FORMAT (///T2s"*DETERMINANT OF TRANSFER MATRIX = "4 
é EAS .a777) 
260 KSUS -="0 
RETURN 
END 
SUBROUTINE EQNS4 (KSUB) 
DIMENSION RMN (25 )sXMN(25) EAR 25925) sNABOR( 2595) 59 
oe AMX(25-.25)3TZNC25) 
CCMMON 7 AIZRMNsXMN SAZ/STZNs AMX /SAS/NASOR 9 EAR 
A JBISAXsDIsHTS /B2/DXsDR /B3B/NXGsNRG /B4/NZsNES93NE 
: /B5/NE1,NES1 
e ZCI/SRMsHTCsESsCPsFLO 


INITIALISE THE MATRIX 
par Fler. NZ 
Cael Jae NZ 

AMR CLPs I) = EARCISs 1) 

2 CONTINUE 
AO = AREA (DIeDI/4) 

ADD THE CONVECTION TERMS (RAD. EQUIVT.) 

(A) GAS ZONES 

Rels "Aes  GAS=SURPFACE 

AW = AO*4*DX/DI 
Delos I =NEST sNzZ 
HCR = HRADQCHTCsTZN(1)) 
M = NABOR(I51)+1 
jaya) J=2 9M 


[IAD ASMA HOME TOT 209 KIATAM ASTEMART* . 287.90) TAMROW SEZ 


< nial acon 
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sr oxt 68 © 
(massarit inesy aeeoee “1st baureanes Rig 


too 
ca: oreo wae ee ma €4) 
(LsDRAs=ay =O 
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wy = Otter? ae 

Qos. 
.* 04 «S408 GSTAJH21 = ROARS Enaas 44°. 0eTN\WND) TARAOS ofr 

(*se *,81 °- 

rETN#TUOTHIRG ABSHD ¢ EMPKS* 2 bey *1*9 TAMROR “pen : 
{NEle** UOTTARATER PF. 

testeas(S?.CrIN) TAMAGS SEL 

(a.F res <asT\) TAMAR ter 
(4,£f50¢.S7T\) TAMROR SEL 

(*: 2vCITUBYRIHOD SAD OADII*—STSAD TAMAOD SEZ 


iN\*224JAD ABRAY . 7 
(\N*e EOTAM SSFQMART FO BSS=VIET* -GEINAND TAMRON TEL 
e* © RIQVAM ASASKART FW TWHAATHEETIO*. BTAA\? YAMADA SE! 
(ANNO, ETS * - 
Oo = BUZy OOS ° 
“wauTan - 
ove 
jauen) eevios ant tyvonave 
o(@e@S ISOGAM. (SS 25 RAS (25 Mes (ERI MMA 4OT2nH3MI0 
(253480 etSS,ASIXMA of 
PAS< SHPAVNZAN SVMASCAETNGAN MME eHMONTAN “MoMM.DD : 
BW e2IW a SUNGEN DAM eDKMNEES ROCKONGGN STs Tt OeKANEBS ~~ « ‘ 
a ee 
Ou eta Ea oD Tee MEY 7 
: > ao 


KIRTAM SHT sargartiMa a. 


: sail Bere ae 
te TVIUGS .GAS) en sat v1 Yoana, SMT a 
: z 7 = ; i. : 


K = NABGR(I.sJ) 
ih GKASGheRNES i) GOTO 5 


AG=TAW 
IF (K LE. NE) A = AREA(2*RMN(K)40R5DR) 
GAL FFSUBS9 (KeIlsFFsi1) 


HC = HCRXFFXxA 
AMX(IsK) = AMX(€1IsK)+HC 
AMXS CPST) = AMXCIs1 )-Hc 


S CONTINUE 
10 CONTINUE 


Cc Aes 


14 AMX(I.51) 


GAS-GAS (BULK FLOW) 
HB = FLO*CP/A0 
DO 20 I=NES1.NZ 
HBR = HRADQ(HBsTZN(I)) 
M = NABOR(I51)+1 
A = AREA(2*RMN(I)4DR,. OR) 
DONT SMJ= 25.M 
K = NABGR(I.J) 
16 €K 2HIGSNES?) GOTO 15 
IF (RMNCK) .NEe RMN(I)) GOTO 15 
IF (XMNCK) «GT. XMNC(I)) GOTO 15 


CALL FFSUBO (IesKsFFs2) 
HBB = HBR*XFFXA 
AMX(1IoK) = AMX(1I4K)+HBB 


AMX(1,1)-H6B 


15 CONTINUE 


Vi= DXx*A 


20 AMX(Io1I) = AMX(I51)-4*ROSK(TZN(1I) 1) %*V 


G Ae4 


mie 
26 


INLET-PLANE GAS ZONES 
I = NES1-NXG 
DO 22 K=1.5.NRG 
I = I+NXG 
A = AREA(2*RMN(1)4DR,s OR) 
AMX€I 51) = AMX(C1IsI) -HRADQ(HBsTZN(I))*A 
+e.) WALL- AND END-— -—ZGNES 
DOe4081=15 NES 
AM=t AW 
TFC GGT. NE) GOTO 26 
A = AREA(2*RMN(I )+DR>,DR) 


AMX(IsI) = AMX(1,1)-CEStHHRADQ(HTSsTZNC(I)))*A 


HCR = HRADQ(HTCsTZN(1) ) 

M = NABOR(I.51)+1 

DOTSS¥ J=2 6M 

KN= NABORCI«J) 

TECK LTS NESTA (GOTGESS 


CALL FFSUBS (IsKeFFs2) 
HC = HCR*FFXA 

AMX CISK)) = AMX€4sK) HC 
AMX(IeI) = AMXCIsTI)-HC 
CONTINUE 

CONTINUE 

KSUBT =" 0 

RETURN 

END 

SUBROUTINE RHSPS (KSUB) 
DIMENSION RHS(25)+RMN( 25) 
COMMON ZAISRMN ZA4/RHS 


/BISAXsDIsHTS /SB2/70XsDR SB3S/NXGsNRG /SB4/NZsNES9 NE 


/B5/NE1 +NES1 


/JCI/RMSHTCsESsCP.FLO /C2/TSURR,TSUPP.CFU,QCOsFLM 


147 


1 
a : 
it 
@i oTos (1234 eta 
#t OTGD (4{1)mMa .34, (aD ne oe 
21 OTOD (illus .Td. Hobo 
(Sy2%e%0t) cave’ a4 
pe en _ 
GEH+ (re TIXMA = (Ned DAMA : 
SeM=(1.TIKMA = pir apse 
BVUMITHOD & 


AeKo eV 
VO¢i se LIMST)AZDACe- Cte TIXMA = CEed os: 
2200S 2A2 JMAINTAUAL A D> | 
ox 125nH © 1) 
oAM.t=4 SS OO. 
oAKeI = 1. 
(8G »Sc+aT JMMROSIASRA = A i . 
At( (145 T+ BHIOGARH—L1,IdRMA = (Ee SOKMA SS. ' 
25403S- =-0M35 OMA -4JAP (3)+fa) > 
234_f=I Oe oa ay 
WA A i 
65 oT0a ( a“ Pe r) ar i 
(OG ,RO+1 DMR" S)ASRA = Aas 
A¥C CCL IYST  STHINAANNG2S1- 1147 I5MA = ted IRMA as 
((2IMST.9THIOGSSH = ADM . 
tel i, 1) A08AM.«- M 
MeS=L 6E GO 
(uel IHOGAM. = Ww 
2E GTQD (£29M ote WP AT. ; 
(Se35. 401) Saue7S JAD 
AFA BH 
Det (Hel DAMA & (hel be 
W-C1 eT XMA = lel 


ya: 


2+ 2aMy SHNOSN 


CoP CAl she END- AND WALL- 
Al AREACRM,RM) 
DO 10 I=1,.NES 
D DI 
DL DX 
LeU Lee Gire tlc t) 
D 2*RMN(I)+4+DR 
DL DR 
S A AREA(D.20L) 
10 RHS(TI) -HTS*A*TSURR 
CORRECT FOR END (EXIT) HOLE PRESENCE 
DCS =e aN 
RI RMN(1) 
IF (RI «GE. RM) GOTO 15 
= A1/AREA(2*RI+DRsDR) 
RHS(I) (1-F)*RHS(I) 
15 CONTINUE 
(Cc) GAS ZONES 
DO 16 I=NES1,NZ 
16 RHS(1) O's 
NOG FLM/DX 
MCO NES4#NOG 
QGN QCOXCFUXDX/FLM 
DO 20 I=NES1.,MCCO 
20 RHS (1) -QGN 
F FLM/DX 
DO 21 J=1.,NGG 
Pants F-1.0 
Lewy CE ~ Glamor ores (MCOs i,) 
INLET-PLANE GAS ZONES 
I = NES1—-NXG 
AO AREAC(DI+-DI/Z4) 
HB = FLO*CP/AO 
DO 23 K=1.,NRG 
lea TENXG 
A = AREA(2*RMN(CI)4DRsDR) 
23 RHS(I) RHS(1)-HB*TSUPP*A 
PRINT 25 
PRINT 305 (RHSC(1),.,T=1sNZ) 
25 FORMAT (///T2s*RHSP5 OUTPUT : 
30 FORMAT (10£13.4) 
CSUR 6 =/.0 
RETURN 
END 
SUBROUTINE SOLV6 (CITR+sERFsKSUB) 
CCMMON 


=ZONES 


SOTO. Ss 


‘/T2 


os RHS OF THE, EQNS 


-F*xQGN 


Uy a | 


SZA2STZNx AMX SA47RHS SASC/SEZN 


1 /B4/NZsNESsNE /SBS/NE1 9 NES1 


DIMENSIGQN 
2 MAJOR( 25) +,KRO( 25) 5 
A AMXV(10510)s 


AMX (25 025) sRHS(25)TZN(25) »EZN( 25) »CFA( 25) > 


= AMT (25 925) sFLU( 25325) »RCF (25) sNPV( 25) 


3 
5 
L=1 NZ 
TE Gl EO, sNELVLOR. 1) .EOQ. NES1)) 
PRINT 2s (AMX(1IsJ)sJ=19NZ) 
14 CONTINUE 


Pein T 
PRINT 
DO 14 


CG 
€ MOVE MATRIX INTO WORKING SPACE 


DBO g20y L=1is:NZ 


(1) 


TEST FOR ILL-CONDITIONING OF THE MATRIX 


PRINT 4 


®AMX* 
‘AMT? 
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A sey pp ina ohn 


Aaa 


MIA ROUSSE ODO = HOO) 
DOM, (2aM=t OS CO 
 4DO- = (1) @4n OS 

xO = 
oOM.t=t ©S oa 
o.1-3 = 7 45) ’ 
Veoss=- = (1e00Mparer {| .0 ato. 4Y Mt? oo oS 
23NGS ZAR SWAUS-TSMI 2 
oxw= Les =". " 
(@¢\IG,yIQIABSA = OA . 
oAN99e0 UF = BH - 
om, f=) ES 00 
oxwet @& Ff.) 7 
(60.901! OMReSTABRA @ A 
AvSGUZTSOM~ (a pee = @2)enR ES — 1) 
@s Turana 
(nM. f=Tel tT see «OE THIRS | © 
(\*? 2403 SHT 320 SHA*sSINeS TUSTUO SH2HA*.STNA\) TARROR ES 
(4.51208) TAMADA OF 


o = 8U2H 
vay 
(Bue A, WIehT LH) SVIG8 2uLTYOReUa 
WSANZA\. SHANRAN 2MAQMETNGAN nOMMID 


Tea. IBANCEN 424-254 Sh \agy 

(ES 74ND « (25 MSAe LOSI AST 4 (ES 12Hity (SS. SS RMA 
» |25 908K, (2598 ad 

+f fe OLIVEMA 

(25 1VS4e t28)- 55h (SSe8S )Uste (Po. 8S 


° . Ubesv 2. 


a 


' ’ be 
2 
- 


tee 


ia ire 
+ 
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DG) Aisy Asal Gin 149 
15 AMT(CIsJ) = AMX(1I.4J) 
20 CONTINUE 
(2) NORMALISE 
BOP s9 St=) NZ 
RMS = 0. 
DOVS5S9=1.NZ 
25 RMS = RMS+AMT( I,J) **2 
RMS = SQRT(RMS/NZ) 
DO: 297 39=15NZ 
29 AMT(IsJ) = AMT(IsJ)/RMS 
30 CONTINUE 
(3) FIND DETERMINANT OF NORMALISED MATRIX 


p= OF 
ME=tG 
GAUL CSLND(AMTs25e9NZsDsMs AMT) 


PRINT 1s Ds™M 
(4) TEST THE MAGNITUDE OF THE DETERMINANT 
i] CABS COs VOe%2M) «GE. 0.5) GOTO 33 
1 SC -CONDIVIGONED USE THE’ *CSLISB SUBROUTINE CSLNIL 
PRINT 4 
PRINT 6 
DOr 2 T=1sNES 
DOWS199=1s5NES 


31 AMXV(IsJ) = AMX(IsJ) 

32 CONTINUE 
GALI CSLNILCAMXVs RHSs EZNs NZs FLUs RCFe NPV) 
GOTO 83 


IF NOT ILL-CONDITIONED TEST FOR DIAGONAL DOMINANCE 
33 MU = 9 
DOS Set — ila NZ 


AI = ABSCAMX(CI,1)) 
HU = TAI 

M =o 

SU = We 


00 §34 =tJ=1.NZ 
if Cf SEG6"E) “GoTo 34 
AJ = ABS(CAMX(I,J)) 
SU = SU+tAJ 
BRECON I Cle Ss Huy YGOTO 34 
HU = AJ 
MoS * J 
34 CONTINUE 
MAJOR (1) = M 
PETCA NEALE SSW) AM = Ms 1 
35 CONTINUE 
TE) OMU SEQ. 0) GOTO 365 
PRINT 4 
PRINT 7 
PRINT 8 
ORINTMA ys SGC TMA IGRGI yet I=1sNZ) 
GOTO 40 
IE SATISFIEDs USE *CSLIB SUBROUTINE CSLNS 
S36. Dees -Oe 
Marana 
GALE CSLNSCAMXs259NZ sRHS 2EZNoDoM AMT ) 
GOTO ''83 
Tf Notes USE THE GAUSS-—SEIDEL ITERATIGN 
(1) COMPUTE THE ORDER OF SOLUTION 
AO HU" = SO) 


XIATAM OBS1IAMHOM AO 1 


(TMA ¢MeGeSMyOSeTMA ION I2> , . A> 
; M20 ny a AT 93 4, 
TAAVIMASTSC ShT WO BOQUTIMOAM SHwT veaT 
EE TOD (8,0 .So. (MeeoTed>es) 7 a 
JIMia2> 241TYOReU2 BIJ2D* SHT Heu OSHOTT a2 — 


2aM,ret 8€ 00 
(LsTIe8A = (usr vKwa KE a 


SUNITMOD SE 
(VQ .43R eUl9 oSH sHSR +2H5 *VKMADITMIEO - Eh no 
AIMAWINOG JSAMODATO #08 TesT manda Preside al ra ae | a 
o = UM Se 
SMsi=t #2€ OO 
(1iel IXMADZBA = TA Pa 
“TA @ UH 
lew 
«0 = ve 
Mei=t BE OC 
aE oTCe (1 .O3. u)  ) in 7 
(ft.l )XMA}BBA = LA - 
Late = Ue 
se oTrd> (um .3y, CAT SI 7” 
LA = UK” 
t=h _ 
SUMITMOD 46 
w = 10) #OLAM ? 


t+ue = UM (Ue ,a4,-TAY ST 
auNTTHOD 

@e c700 (6 03, a a 
THIAG 

: THER 


USVel=1 of fPROLAMe TDS ot Bate 
96 zo 


24323 SHITVORSUe Siljede sau a1 Ae 


he 


DOmSo late NZ 
AJ = ABSCAMX(NZsJ)) 
if GAdy .Ue. HU) «GOTO 38 
HUy = AJ 
M=J 
38 CONTINUE 
MAJGR(NZ) = M 
has «NZ 
0D1 60 pwatil=2.NZ 
PAs LA 
I =, t=1 
B04 S. K=1,NZ 
45 KROGK)« = ol 
K = NZ+1 
DO 450° KK=1e TA 
K = K-1 
50 KRO(CMAJOR(K)) = O 
HU A= «0.5 
BOSS J=1.NZ 
AJ = KRO(J)*ABS(CAMX(IsJ) ) 
ihe AS ecole chu). .GOTO 59 
TU AS 
M =<e«J 
59 CONTINUE 
60 MAJOR(I) = M 
DOGG FI =1eNZ 
61 NPVCI) = MAJOR(NZ-141 ) 
PRINT 9 
PRINT SlOs,GNEVGhis L=teNZ) 
(2) SOLVE EQNS IN THE PRE=-DETERMINED ORDER 
bGpGLarR peGleublm GOTOW.66 
DGeacsS L=psNZ 
65 EZNGLp =/CGNVRTCEZN(I) s1) 
66 1 = NZ+1 
DO 80 I1I=1»sNZ 
LAsel=1 
SOunOekK=1.NZ 
70 KRO(K) = 1 
GQn75 K=isI 
75 KROCMAYJOR(K)) = O 
SU = RHS(I) 
Gein Sar s= 15 NZ 
79 SU = SU-KRO(J)*AMX( 19d) *EZN(J) 
89 EZNCMAJOR(I)) = SU/AMXCIsMAJOR(I)) 
COMPUTE NEW TEMPS AND ERROR PARAMETER 
83 Che= F150 
MUic= 4 
POVSSn1=15.NZ 
Hem CEeZNGl) Le. ivis. -)) “MU = MUTI 
85 CONTINUE 
bey CMU) se0s6: Ady GOTOed 21 
Milo d 
DOpi1 0, i=l NZ 
EO = CONVRT(TZN(I)91) 
110 CFA(1) = ABS(EZNCI1)/EQ-1.0) 
GEusy CEA 1) 
DG = 11S 1 =25.NZ 
Cle 9GrnGt) 
LEaCedes bowed «em GDite ls 
Peel boieeCh.feCreeeat 
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elitr 1) 
I=1,NZ 


Lele 0.5*CF 


CONVRT(TZN(1)51) 
CF*EZN(1)+(€1-CF) *E0 


gO. Lri4+. ) GOTO 1206 


MU+1 
21 


CONT INUE 
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T=15,NZ 


CF*xEZN(1)4+€1-CF) *CONVRTCTZNCI) 91) 
CONVRT CEN 92) 
ERF+(T-TZN(1) ) *x2 
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CONTINUE 


PRINT 1359 
PRINT 12% 
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SUBROUTINE 
ENERGY BALANCE FOR ENTIRE SYSTEM 


COMMON 


° /BisAXs,DI 
° /B4/NZ NES eNEsNSs NG 

; 7CI/SRMsHTCsESsCPsFLO 
DIMENSIGN 


1 


MU 
CF 


SQRT (ERF/NZ) 


FORMAT (///T2s*DETERMINANT OF NORMALISED MATRIX =", 
. F10.5,°E*,13) 
FORMAT (/10E13.4) 
FORMAT (*1°*) 
FORMAT (//1X) 
FORMAT (//T40s"*COEF MATRIX OF THE ENERGY BALANCE EQNS#//) 
FORMAT (//T2s*COGEF MATRIX ILL-CONDITIGNED USE CSUNTi) 
FORMAT (//T25s"COEF MATRIX NOT OLAG-DOMINANT*/T2s5 
a "PROCEED WITH GAUSS-S=EIDEL ITERATION!®/) 
FORMAT (/T2s*THE DOMINANT TERMS ARE 3 '*/) 
FORMAT (//T2s*THE ORDER OF SOLUTION IS :*7) 
FORMAT (¥T20>3 515) 
FORMAT (/7T255(21493X)) 
FORMAT (//T2s "MAX CORRECTION FACTOR = *,£11.4/) 
FORMAT (///T20"°# OF TRIAL CF-VALUES =',110) 
KSUB = 0 
RETURN 
END 
SUBROUTINE QFLX7 (KSUB) 
DIMENSION TZN C25) 5 OF XCeS) 
COMMON JA2/STZN SAB/QFX 
1 /B1I/SHTS /B4/NZsNES 
2 /C1/RMs HTC /C2/TSURR 
DOO ei—15 NES 
QFX(I) = HTS*(TZNCI)-TSURR) 
KSUB = 0 
RETURN 
END 
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AS “A CHECK ON CALCS 

SA2/TZN SAB/QFX /SAG/EZN 
/B3/NXGsNRG 

/BS/NE1NES1 
/C2/TSURR»s TSUPP 3CFUsQCOsFLM 

TZN(25) > QFX(25) sRMN(25) » 

EZN(25) 


ZALISRMN 
7B2/DXsDR 


NET ENERGY GENERATION RATE 


QGEN 


QCO*CFU-CP*FLO*(TZNCNES#NXG) -TSUPP) 
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C TOTAL NET HEAT FLUX (SURFACE ZONES) 
GEEUX = 0. 
DO 5 I=1sNEs2 
A = AREA(2*RMN(1I)40R,0R) 
5 QFLUX = QFLUX+A*QFX(I) 
Al = AREA(RMsRM) 
OOTiGw T=2,NE.2 
RI = RMN(T) 
A = AREAC2*RI+DR.OR) 
TRVCRI #SERS RM) CLAC=EPA—A1 
QFLUX = QFLUX+A*QFX(I) 
190 CONTINUE 
A = AREA(DI.0x) 
DO 15 I=NE1.,NES 
1S QFLUX = QFLUX+A*QFX(TI) 
C COMPARE 
PC = QFLUX/QGEN*100. 
C OUTPUT 
PRINT 50s QFLUXs PCs QGEN 
50 FORMAT (C7/777TZ2ZSCtEBLGCSsQUTPUT 2*/7/T2s"TOTAL NET FLUX =", 
eo “EVSS6s5X"*BTU/SHR*//T2s "WHIGH 15*sF8.202X%09'% GE*//T2s 
e “THE NET ENERGY GEN RATE OF '5E15.655Xs"BTU/SHR® ) 


KSUB = 0 

RETURN 

END 

SUBROUTINE FFSUB9O (IsJ5FFsK) 
DIMENSION TZN(C25) 

COMMON SA2/TZN 

Til = TZNERIAT ZN) 

Sas TiueTt 

FF = 1/((14+T1)*C€14+T2) ) 

PRee CK SE0nee) FE = T2*T1*FF 
RETURN 

END 


FUNCTION CONVRT(Qs5K) 
C CONVERTS BETWEEN TEMP AND EMISSIVE PWR 
C --DEPENDING ON K 

COB=QT7T4 50 

HG= L000 0 


Po=14.6 
TEBCKSCEG. 1) GOTO,” S 
HGo= co 
co = D 
Dp = H 
Pas. 0.25 
S CONVRT = CO*(Q/D)**P 
RETURN 
END 


FUNCTION AKGR(TEMP) 

C COMPUTES SUM OF WT-ING FACTOR-ABSORPTION COEF PRODUCTS 
DIMENSION B8N(5510)5SK(5) 
COMMGN /ZAS3/BNsSK /C2/NGRAYsNB 
AKGR#=804 
PO L0ATE1 sNGRAY 
HU? 150 
AG = BN(Is1) 
ony, Sey=2.NS8 
T = TxXTEMP 

5S AG = AGtBN(I1I5J)*T 
10 AKGR = AKGR+AG*SK(I) 
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FUNCTION HRADQ (HoT) 
C COMPUTES AT TEMP T PRODUCT OF AREA AND RADIATIVE =ZQUIVT 
C OF CONV HT TRANSFER COEF 
HRADQ = H*XT/SCONVRT(T)1) 
RETURN 
END 
FUNCTION AREA (OIA,DEL) 
C COMPUTES AREAS OF CIRCLES, CIRCULAR ANNULI AND CYLINDERS 
PI = 0.31415926E 01 
AREA = PI*xDELXDIA 
RETURN 
END 
FUNCTION ROSK (TEMPsKRO) 
GIVES THE POSSELAND WEAN ABSQRPTIGN COEF 
GR A RELATED MEAN ABSORPTION COEF 
(DEPENDING ON KRO) 
FOR THE GRAY GAS MIXTURE 
DIMENSION 8N(5519)+SK(5) 
COMMON /A2/BNsSK /C2/NGRAYs NB 
ROSK = =..0.< 
DO 10 I=1,s,NGRAY 
Tae le 
SU = BN(I.o1) 
DO 5S J=2,NB 
T = T*TEMP 
p= J+s 
EJP On25%eJ 
5 SU = SUtFI*BN(I 55) *T 
10 ROSK = ROSK+SU/(SK(I1) **KRC) 
ROSK = 1/ROSK 
RETURN 
END 
FUNCTION AGO (TEMP) 
C CALCULATES AT TEMP T CLEAR GAS WT-ING FAC 
DIMENSTON BN(5510) 
COMMON /A3/BN /C3/NGRAY>sNB 


(a) a) @ 


AGO = 1.0 
DO 10 TI=1s,NGRAY 
T = 1.0 


AGO = AGO-BN(I.1) 
DO 5S J=2.NB 
T = T*TEMP 


5 AGO = AGO-BN( IJ) *T 
10 CONTINUE 
RETURN 
END 
FUNCTION DIRECT(AsB2C) 


G GOMPUTES DIRECT EXCH AREAS BETWEEN 
Cc ANY PAIR OF "GCOAXKAL CIRCULAR SUREAGES 
Pi w= O.oLouscoeZe EO 
ABC = A*A4+B*B4C#*C 
AB = Ax*B 
SUBT = SQRT(AECXABC—4 *XAB*¥AB) 
DIRECT = 0.5*PI*( ABC-SUBT) 
RETURN 
END 
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